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ABSTRACT 

We present new limits on an isotropic stochastic gravitational-wave background (GWB) using 
a six pulsar dataset spanning 18 yr of observations from the 2015 European Pulsar Timing 
Array data release. Performing a Bayesian analysis, we fit simultaneously for the intrinsic 
noise parameters for each pulsar, along with common correlated signals including clock, and 
Solar System ephemeris errors, obtaining a robust 95% upper limit on the dimensionless strain 
amplitude A of the background of A < 3.0 x 10“^^ at a reference frequency of lyr“^ and a 
spectral index of 13/3, corresponding to a background from inspiralling super-massive black 
hole binaries, constraining the GW energy density to Q.^^{f)h^ < 1.1 x 10“^ at 2.8 nHz. We 
also present limits on the correlated power spectrum at a series of discrete frequencies, and 
show that our sensitivity to a fiducial isotropic GWB is highest at a frequency of 5x10“^ Hz. 
Finally we discuss the implications of our analysis for the astrophysics of supermassive black 
hole binaries, and present 95% upper limits on the string tension, G///c^, characterising a 
background produced by a cosmic string network for a set of possible scenarios, and for a 
stochastic relic GWB. For a Nambu-Goto field theory cosmic string network, we set a limit 
Gfilc^ < 1.3 X 10“^, identical to that set by the Planck Collaboration, when combining Planck 
and high-^ Cosmic Microwave Background data from other experiments. For a stochastic relic 
background we set a limit of < 1.2 x 10“^, a factor of 9 improvement over the most 

stringent limits previously set by a pulsar timing array. 
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1 INTRODUCTION 


The first evidence for gravitational-waves (GWs) was originally ob¬ 
tained through the timing of the binary pulsar B1913-I-16. The ob¬ 
served decrease in the orbital period of this system was found to 
be completely consistent with that predicted by general relativity, 
if the energy loss was due solely to the emission of gravitational 
radiation jTaylor & Weisberg||1989^ . Despite a decrease of only 
2.3ms over the course of 30 yr, by exploiting the high precision 
with which the time of arrival (TOA) of electromagnetic radiation 
from pulsars can be measured, deviations from general relativity 
have been constrained by this system to be less than 0.3% jWeis-| 
|berg, Nice & Taylor|2010| l. 

Since then, observations of the double-pulsar, PSR 
J0737-3039, have provided even greater constraints, placing 
limits on deviations from general relativity of less than 0.05% 
( [Kramer et al.| j2006| l, Kramer et al. in prep.). It is this extraor¬ 
dinary precision that also makes pulsar timing one possible route 
towards the direct detection of GWs, which remains a key goal in 
experimental astrophysics. 

For a detailed review of pulsar timing we refer to jLorimer &| 
[Kramerj ( |2005[ (. In general, one computes the difference between 
the expected arrival time of a pulse, given by a pulsar’s timing 
model which characterises the properties of the pulsar’s orbital mo¬ 
tion, as well as its timing properties such as its spin frequency, and 
the actual arrival time. The residuals from this fit then carry phys¬ 
ical information about the unmodelled effects in the pulse propa¬ 
gation, including those due to GWs (e.g. |Sazhin|197 8l|Detweiler| 
[TW9l l. 

Individual pulsars have, for several decades, been used to set 
limits on the amplitude of gravitational radiation from a range of 
sources (e.g. [Kaspi, Taylor & Ryba|[T994) . However, by using a 
collection of millisecond pulsars, known as a pulsar timing array 
(PTA, [Foster & Backer|ri990| l, one can both increase the signal- 
to-noise ratio of the etfect of gravitational radiation in the timing 
residuals, and use the expected form for the cross correlation of the 
signal between pulsars in the array to discriminate between the GW 
signal of interest, and other sources of noise in the data, such as the 
intrinsic spin-noise due to rotational irregularities (e.g. [Shannon &| 
[C ordes[2010[ l, or delays in the pulse arrival time due to propagation 
through the interstellar medium (e.g. [Keith et al.|2013^ . In the spe¬ 
cific case of an isotropic stochastic gravitational-wave background 
(GWB), which is the focus of this paper, this correlation is known 
as the ‘Hellings-Downs’ curve ( [Hellings & Dow^[1983) , and is 
only a function of the angular separation of pairs of pulsars in the 
array. 


The lowest frequency to which a particular pulsar timing 
dataset will be sensitive is set by the total observing span for that 
dataset. Sensitivity to frequencies lower than this is significantly 
decreased due to the necessity of fitting a quadratic function in the 
pulsar timing model describing its spin down. PTA datasets are now 
entering the regime where observations span decades, and as such 
are most sensitive to GWs in the range 10“® - 10“® Hz. The primary 
GW sources in this band are thought to be supermassive black hole 
binaries (SMBHBs) ([Rajagopal & Romani[[1995 JJaffe^^^acto 


2003|[Wyithe & Loeb|2Q03[ [Sesana et al.|2004[ 


Sesana, Vecchio 


& Colacino| 2008TJiowever other sources such as cosmic strings 

(see, e.g. [VilenkinJ 1981 [[Vilenkin & Shellard|1994[ l or relics from 
inflation (see, e.g. lGrishchuk 20051 have also been suggested. 

The formation of SMBHBs is a direct consequence of the hi¬ 
erarchical structure formation paradigm. There is strong evidence 
that SMBHs are common in the nuclei of nearby galaxies (see[Kor-| 


[mendy & Ho[2013| and references therein). The fact that many dis¬ 
tant galaxies harbour active nuclei for a short period of their life im¬ 
plies that they were also common in the past. In A-Cold Dark Mat¬ 
ter (A-CDM) cosmology models galaxies merge frequently (|Lacey| 
[& Cole|1993^ . During a galaxy merger the SMBHs harboured in the 
galactic nuclei will sink to the center of the merger remnant, even¬ 
tually forming a SMBHB ( [Begelman, Blandford & Rees 1980) . As 
a consequence the Universe should contain a potentially large num¬ 
ber of gradually in-spiralling SMBHBs. The incoherent superposi¬ 
tion of GWs from these binaries is expected to form an isotropic 
stochastic GWB. Deviations from isotropy, however, such as from 
a small number of bright nearby sources, could result in individ- 


ually resolvable systems (IL 

ee et al.|2011 1, and an anisotropic dis- 

tribution of power across the sky (Mmgarelli et al.|2013| [Taylor &[ 

[Gair|2013|[Gairetal.|2014 

1 . These latter situations are the subject 

of two companion papers ( 

Taylor et al.[[2015[ [Babak et al.|2015[l 


here we focus on the possibility of detecting a stochastic isotropic 
GWB, and we will discuss the implications of our findings for the 
astrophysics of SMBHBs, cosmic strings, and relics from inflation. 


An isotropic, stochastic GWB of cosmological or astrophysi- 
cal origin can be described in terms of its GW energy density con¬ 
tent Pg„ per unit logarithmic frequency, divided by the critical en¬ 
ergy density, pc, to close the Universe: 


Dgw(/) = 


1 dpg,v 

Pc din/ 


27r^ , , 


( 1 ) 


Here, / is the GW frequency, pc = is the critical energy 

density required to close the Universe, Hq = 100 h km s“* Mpc“* is 
the Hubble expansion rate, with h the dimensionless Hubble param¬ 
eter, and pg„ is the total energy density in GWs ([Allen & Romano [ 
[19991 [Maggiore[2000[ l. 

Typically the ‘characteristic strain’, hdf), associated with a 
GWB energy density Dg„(/) is parametrised as a single power-law 
for several backgrounds of interest: 



( 2 ) 


where A is the strain amplitude at a characteristic frequency of 
lyr^', and a describes the slope of the spectrum. Finally, he is di¬ 
rectly related to the observable quantity induced by a GWB in our 
timing residuals, the one-sided power spectral density, S (/), given 
by: 


S(f}- 


1 1 


1 

(^ y 

12;r2 ' 

Ur-'j 




yr 


(3) 


where y = 3 - 2a. Note that unless explicitly stated otherwise, 
henceforth when referring to spectral indices we will be referring 
to the quantity y. 

The expected spectral index varies depending on the source of 
the stochastic background. For a GWB resulting from inspiraling 
SMBHBs the characteristic strain is approximately heif) oc 
([Rajagopal & Romani|1995[[Jaffe & Backer[2003[[Wyithe & Loeb[ 
[2003[ [Sesana et al.|2004^ ,' or equivalently, y = 13/3, whereas pri¬ 

mordial background contributions or cosmic strings are expected to 
have power-law indices of y = 5 ([G rishchuk|2005[l, and y = 1 6/3 

( Olmez, Mandic & Siemens|2010 Damour & Vilenkin||2005| re¬ 

spectively. However, for cosmic strings in particular, a single spec¬ 
tral index is not expected to accurately describe the spectrum in the 
PTA frequency band |Sanidas, Battye & Stappers[2012^ . 

A multitude of experiments have set limits on the amplitude 
of the stochastic GWB, either at a reference frequency as is done 
for PTAs ([Shannon et al.[[2013[l and ground-based interferome- 
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Figure 1. Summary of key results from the analysis of a 6 pulsar dataset from the 2015 EPTA data release (D15). Results are presented in terms of Qg^C/) as a 
function of GW frequency, with Hq = 70km Mpc“^. We indicate the 95% upper limits on the amplitude of a con'elated GWB assuming a power law model 
with a spectral index of y = 13/3 (solid black line; Section[5j and for a more general analysis where the power is determined simultaneously at a set of discrete 
frequencies (dashed line) as discussed in Section [5.1.1| The red shaded areas represent the central 68%, 95%, and 99.7% confidence interval of the predicted 
GWB amplitude according to |Sesan^ l |2013b) under the assumptions that a SMBHB evolves purely due to gravitational radiation reaction and binaries are 
circular (See Section[6^for more details). Only about 5% of the distribution is excluded, meaning that our limit does not place significant restrictions on the 
cosmic SMBHB population. We also indicate 95% upper limits obtained for a stochastic relic background (green dash-dotted line; Section [63) , and for cosmic 
string network backgrounds (blue triple-dashed line; Section [6!^ . The cosmic string limit plotted corresponds to a fiducial model for a population of cosmic 
strings, with the following parameters: string tension Gfijc^ = 10“^, the birth-scale of loops relative to the horizon o^cs = 1-6 x 10“^, spectral index q = 4/3, 
cut-off on the number of emission harmonics n* = 1 and intercommutation probability /> = 1. Finally we indicate recent constraints placed by CMB jSendra] 
|& Smith|2012) , and BBN jAllen|1997||Maggiore|200^|Planck Collaboration et al.|2015) observations. 


ters jAasi et al.|2014| l, or by reporting a value for GW energy den¬ 
sity integrated over all frequencies as is done by Big Bang Nu- 


cleosynthesis measurements, e.g. ^Cyburt et al. 

2005} and Cosmic 

Microwave Background (CMB) measurements 

Smith, Pierpaoli &| 


|Kamionkowski|[T006[ |Sendra & Smith|[2012^ . As such, an upper 
limit on the stochastic GWB reported in terms of either 
or flgwl/) for a specified value of h provides a clear way to report 
our limits. 


In the last few years the European PTA (EPTA), Parkes PTA 
(PPTA), and the North American NanoHertz Observatory for Grav¬ 
itational waves (NANOGrav) have placed 95% upper limits on 
the amplitude of a stochastic GWB at a reference frequency of 
lyr“* of 6x10“'^ jvan Haasteren et al. 2011 1 , 2.4x10“*^ jSh; 


I et al.pofs ', and 7x10 (Demorest et al.|2013l respectively. 


While many of the same pulsars are used by all the PTAs, and 
both the EPTA and PPTA have similar total observing spans, all 
of these limits have been placed using different datasets, and dif¬ 
ferent methodologies. As such, these similarly constraining limits 
should not be seen as redundant, but rather as complementary. Eor 
example, the first EPTA limit used Bayesian analysis methods, pro¬ 
ducing an upper limit while simultaneously fitting for the intrinsic 
timing noise of the pulsars. Subsequent limits have used simula¬ 
tions to obtain conservative upper bounds consistent with the data. 


or made use of frequentist methods, fixing the noise at values de¬ 
rived from analysis of the individual pulsars. Naturally a simultane¬ 
ous analysis of the intrinsic properties of the pulsars with the GWB 
is the preferred method, and we will show explicitly in Section]^ 
that fixing the noise properties of the individual pulsars can lead to 
an erroneously stringent limit on the amplitude of a GWB in the 
pulsar timing data. The three PTA projects also work together as 
the International PTA (IPTA; [Hobbs et al.||2010^ , where all three 
datasets are combined in order to produce ever more robust and 
constraining limits on the GWB, with the eventual goal of making 
a first detection. 

In this work we make use of the Bayesian methods presented 
in jLentati et al.| ( |2()13^ (henceforth L13), which allows us to greatly 
extend what is computationally feasible for a Bayesian analysis of 
pulsar timing data. In particular, while we obtain upper limits on 
the amplitude of a GWB using a simple, two parameter power law 
as in jvan Haasteren et al.H20lT} , we can also make use of a much 
more general model, enabling us to place robust limits on the corre¬ 
lated power spectrum at discrete frequencies. We can also include 
additional sources of common noise in our analysis simultaneously 
with the GWB, such as those that could be expected from errors 
in the Solar System ephemeris, or in the reference time standard 
used to measure the TOAs of the pulses. Einally we also take two 
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approaches to parameterising the spatial correlations between pul¬ 
sars, without having to assume anything about the form it might 
take. This spatial correlation is the ‘smoking gun’ of a signal from 
a GWB, and so the ability to extract it directly from the data is cru¬ 
cial for the credibility of any future detections from pulsar timing 
data. 


The key results of our analysis, compared to current theoreti¬ 
cal predictions for a range of models of stochastic background and 
indirect limits in the PTA range are summarised in Fig. ^ In Sec¬ 
tion we describe the deterministic and stochastic models that we 
included in this analysis. In Sectionj^we discuss the implementa¬ 
tion of these methods in our Bayesian and frequentist frameworks. 
In Section|^we introduce the EPTA dataset adopted for the analysis 
(Desvignes et al. in prep.), and in Section we present the results 
obtained from our analysis. The implications of our findings for 
SMBHB astrophysics, cosmic strings, and relics from inflation are 
discussed in Section]^ and finally, in Section|^we summarise and 
discuss future prospects. 

This research is the result of the common effort to directly 
detect gravitational-waves using pulsar timing, known as the EPTA 
(EPTA; [Kramer & Champion|2013| in 


2 SIGNAL AND NOISE MODELS 

The search for a stochastic GWB in pulsar timing data requires 
the estimation of a correlated signal of common origin in the pulse 
TOAs recorded for the different pulsars in the array. The difficulty 
lies in the intrinsic weakness of the signal and the presence of a 
range of effects - both deterministic and stochastic - that conspire 
to mask the signal of interest. At the heart of our analysis methods 
is the variance-covariance matrix 

'v,Ai,j] = {djmAn), ( 4 ) 

that describes the expectation value of the correlation between TOA 
i from pulsar /, with a TOA j from pulsar J. In the following 
description upper-case latin indices identify pulsars, and 

lower case latin indices i, j,... are short hand notation for the TOAs 
t„t„....Eq.g depends on the unknown parameters that describe 
the model adopted to describe the data and enter the likelihood 
function in the Bayesian analysis, and the optimal statistic in the 
frequentist approach. 

For any pulsar we adopt a model for the observed pulse TOAs, 
which we denote d, that results from a number of contributions and 
physical effects according to: 


• T*''’*, the stochastic contribution due to ‘common noise’, 
present across all pulsars in the timing array (described in Sec. |2.6[ l, 
as could be expected from errors in the Solar System ephemeris, or 
in the reference time standard used to measure the TOAs of the 
pulses. 

• T*'"', the stochastic contribution due to a GWB. 

Our model assumes that all stochastic contributions are zero 
mean random Gaussian processes. Each of the contributions just 
described depends on a number of unknown parameters that need 
to be simultaneously estimated in the analysis. While all these el¬ 
ements, which we set out in detail below, will be present in the 
Bayesian analysis described in Section [3T| we do not incorporate 
the common pulsar noise terms in the frequentist optimal-statistic 
analysis described in Section as this approach by design in¬ 
terprets all cross-correlated power as originating from a stochastic 
GWB. 

2.1 The timing model 

The first contribution to the total signal model that we must con¬ 
sider is the deterministic effect due to the intrinsic evolution of 
the ‘pulsar clock’, encapsulated by the pulsar’s timing ephemeris. 
We identify with e; the m-dimensional parameter vector for pul¬ 
sar I that describes the relevant set of timing model parame¬ 
ters, and denote as T(e) the set of arrival times determined by 
the adopted model and specific value of the parameters. We use 
Tempo2 ( |Hobbs, Edw ards & Ma nchester||200^ [Edwards, Hobbs| 
|& Manchester|[2006l l to construct a weighted least-squares fit, in 
which the stochastic contributions have been determined from a 
Bayesian analysis of the individual pulsars using the TempoNest 
plugin ( [Lentati et al.|2014) . We can define the set of ‘post-fit’ resid¬ 
uals that results from subtracting the predicted TOA for each pulse 
at the Solar System Barycenter from our observed TOAs as: 

dposi = d - T(e). (6) 

In everything that follows, rather than use the full non-linear tim¬ 
ing model we consider an initial estimate of the m timing model 
parameters eo, and construct a linear approximation to that model 
such that any deviations from those initial estimates are encapsu¬ 
lated using the m parameters Se such that: 

de,- = - £0,- (7) 

Therefore, we can express the change in the post-fit residuals that 
results from the deviation in the timing model parameters 6e as: 

dt = dpost - MSe, (8) 


d = r™ + rWN+^SN + ^M^^CN+^GW^ (5) 

In Eq. 0 we have: 

• T™, the deterministic model that characterises the pulsar’s as¬ 
trometric properties, such as position and proper motion, as well as 
its timing properties, such as spin period, and additional orbital pa¬ 
rameters if the pulsar is in a binary. 

• T™, the stochastic contribution due to the combination of in¬ 
strumental thermal noise, and intrinsic pulsar white noise. 

• the stochastic contribution due to red spin-noise. 

• the stochastic contribution due to changes in the disper¬ 
sion of radio pulses traveling through the interstellar medium. 


www.epta.eu.org/ 


where M is the N^xm ‘design matrix’ which describes the depen¬ 
dence of the timing residuals on the model parameters. 

When we perform our Bayesian GWB analysis, we will 
marginalise analytically over the linear timing model, as described 
in Section lrT] When performing this marginalisation the matrix M 
is numerically unstable. To remedy this issue we follow the same 
process as in jvan Haasteren & Vallisn^ ( |2014| and take the SVD 
of M, to form the set of matrices USV^ . Here U is an NjxNd matrix, 
which we can divide into two components: 

U = (G^,G), (9) 


where G is a Aj x (Nd — m) matrix, which can be thought of as a 
projection matrix (van Haasteren & Levin 2013 i, and G*" is the N^x 
m complement. G'^ represents a set of orthonormal basis vectors 
that contain the same information as M but is stable numerically. 
We therefore replace M with G'" in the subsequent analysis. 
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2.2 White noise 

We next consider the contribution to the total signal model that 
results from a stochastic white noise component, r™. This noise 
component is usually divided into two components, and this is the 
model that we adopt in our analysis: 

• For a given pulsar /, each TOA has an associated error bar, 
<T(/,,), the size of which will vary across a set of observations. We 
can introduce an extra free parameter, referred to as EFAC, to 
account for possible mis-calibration of this radiometer noise. The 
EFAC parameter therefore acts as a multiplier for all the TOA error 
bars for a given pulsar, observed with a particular ‘system’ (i.e. a 
unique combination of telescope, recording system and receiver). 

• A second white noise component is also used to represent 
some additional source of time independent noise, which we call 
EQUAD, and adds in quadrature to the TOA error bar. In princi¬ 
ple this parameter represents something physical about the pulsar, 
for example, contributions from the high frequency tail of the pul¬ 
sar’s red spin-noise power spectrum, or jitter noise that results from 
the time averaging of a finite number of single pulses to form each 
TOA (see e.g. [Cor des & Downs||1 985[ |Liu et af]|20 11[ | Shannon] 
|et al.|201^ . While this term should be independent of the observing 
system used to generate a given TOA, differences in the integration 
times between TOAs for different observing epochs can muddy this 
physical interpretation. 

We can therefore modify the uncertainty cry ,), defining <T(/_,) such 
that the statistical description is; 

= SuSijdi,,^ ( 10 ) 

where 

^l.i) = («(/,.)0-(/,o)^ +fih) (11) 

where a and represent the EFAC and EQUAD parameters applied 
to TOA i for pulsar I respectively. In Section|^we list the number 
of different observing systems per pulsar used in the analysis pre¬ 
sented in this paper. 


2.3 Spin-noise 

Individual pulsars are known to sometimes suffer from ‘spin-noise’, 
which is observed in the pulsar’s residuals as a red noise process. 
This is a particularly important noise source, as most models for 
a stochastic GWB predict that this too will induce a red spectram 
signal in the timing residuals. The spin-noise component is specific 
to each individual pulsar, and is uncorrelated between pulsars in the 
timing array. The statistical properties of the spin-noise signal are 
therefore given by: 

(rfMTf[;]> = d„C™y), (12) 

where the matrix element C®/’* denotes the covariance in the spin- 
noise signal between residuals i, j for pulsar /. In order to construct 
the matrix we will use the time-frequency method described 
in LI3, which we will summarise below. 

We begin by writing the spin-noise component of the stochas¬ 
tic signal as 

= F^'^asN (13) 

where the matrix denotes the Fourier transform such that for 
signal frequency v and time t we will have both: 
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F^^(v, l) = sin (27rvr), 


(14) 


and an equivalent cosine term, and asN are the set of free parameters 
that describe the amplitude of the sine and cosine components at 
each frequency. 

We include in our model the set of frequencies with values 
n/r, where T the longest period to be included in the model and 
the number of frequencies to be sampled is nsN- In our analysis 
presented in Section]^ we take T to be ~ 18 yr, which is the to¬ 
tal observing span across all the pulsars in our dataset, and we 
take nsN = 50 such that we include in our model periods up to 
~ 130 days which is sufficient to describe the stochastic signals 
present in the data (Caballero et al. in prep.). For typical PTA data 
|Lee et al!]p012| l and |van Haasteren & Levin|P013| l showed that 
taking T to be the longest time baseline in the dataset is sufficient 
to accurately describe the expected long-term variations present in 
the data, as the quadratic term present in the timing model signifi¬ 
cantly diminishes our sensitivities to periods longer than this in the 
data. 

The covariance matrix of the spin-noise coefficients a*'’’ be¬ 
tween pulsars I, J at model frequencies i, j, which we denote 
will be diagonal, with components: 


'P 


SN 

(l.J.'.j) 


/„SN „SN \ _ ,„SN r r 

~ Vl,i 


(15) 


where the set of coefficients represent the theoretical power 
spectrum of the spin-noise signal present in pulsar I. In our analysis 
of the dataset presented in Section]^ we assume that this intrinsic 
spin-noise can be well described by a 2-parameter power law model 
in frequency, given by: 


,^"”(y,AsN,rsN)= 



[ ' ] 

^ 3 y-rsN 

12;r2 1 

ilyrj 

1 y 


(16) 


with Asn and ysN the amplitude and spectral index of the power 
law. 


We note that as discussed in LI3, whilst Eq. {15} states that 
the spin-noise model components are orthogonal to one another, 
this does not mean that we assume they are orthogonal in the time 
domain where they are sampled, and it can be shown that this non- 
orthogonality is accounted for within the likelihood ([van Haasteren| 


& Vallisneri 


2015 I. The covariance matrix C™ for the red noise 


signal present in the data alone can then be written: 


C™ = N7‘-N7'F™[(F®VN7‘F™-i-('I'™)-‘] '(F™)^N7‘, (17) 


with N/ the diagonal matrix containing the TOA uncertainties, such 
that = d-2 . 


2.4 Dispersion measure variations 

The plasma located in the interstellar medium (ISM) can result in 
delays in the propagation of the pulse signal between the pulsar and 
the observatory. Variations in the column-density of this plasma 
along the line of sight to the pulsar can appear as a red noise signal 
in the timing residuals. 

Unlike other red noise signals however, the severity of the ob¬ 
served dispersion measure (DM) variations is dependent upon the 
observing frequency, and as such we can use this additional infor¬ 
mation to isolate the component of the red noise that results from 
this effect. 

In particular, the group delay tg{Vo) at an observed frequency 
Vo is given by the relation: 

tgivo) = DMKKvl) (18) 
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where the dispersion constant K is defined to be: 

K = 2.41 X 10^‘® pc s^‘ (19) 


and the DM is defined as the integral of the electron density from 

the Earth to the pulsar: 


DM = 


f 


nAl. 


( 20 ) 


While many different approaches to performing DM correc¬ 
tion exist (e.g. |Lee et al.| ( |20T4t ; |Keith et al.| ( |2013^ ), in our analysis 
we use the methods described in LI3. DM corrections can then be 
included in the analysis as an additional set of stochastic parameters 
in a similar manner to the intrinsic spin-noise. Further details on 
the DM variations present in the EPTA dataset, including compar¬ 
isons between different models, will be presented in a seperate pa¬ 
per (Janssen et al. in prep.). In our analysis, as for the spin-noise, we 
assume a 2-parameter power law model, with an equivalent form to 
Eq. 1 16 ', however we omit the factor 127r^ for the DM variations. 

We first define a vector D of length Nd for a given pulsar as: 

A = l/(^<„) (21) 


for observation i with observing frequency 

We then make a change to Eq. l |14^ such that our DM Fourier 
modes are described by: 


F°'^(v, f,) = sin (Invtj) Di 


( 22 ) 


and an equivalent cosine term, where the set of frequencies to be 
included is defined in the same way as for the red spin-noise, such 
that we choose the number of frequencies, «dm, to also be 50. Un¬ 
like when modelling the red spin-noise, where the quadratic terms 
in the timing model that accounts for pulsar spin-down acts as a 
proxy to the low frequency (v < l/L) fluctuations in our data, we 
are still sensitive to the low frequency power in the DM signal. As 
such these terms must be accounted for either by explicitly includ¬ 
ing these low frequencies in the model, or by including a quadratic 
term in DM to act as a proxy, defined as: 

Q°^{ti) = qot,D,+qi^D, (23) 

with ^0 1 free parameters to be fit for, and t, the barycentric arrival 
time for TOA i. This can be achieved most simply by adding the 
timing model parameters DM\ and DM2 into the pulsar timing 
model, which are equivalent to qo and qi in Eq. l |23[ >, and this is the 
approach we take in our analysis here. 

As for the spin-noise component we can then write down the 
time domain signal for our DM variations as: 

= F°“aDM, (24) 


with aDM the set of free parameters that describe the amplitude of 
the sine and cosine components at each frequency. 

The covariance matrix of the coefficients a°^ between pulsars 
1, J at model frequencies i, j, which we denote is then equiv¬ 
alent to the spin-noise matrix in Eq. {ID, and we can similarity 
construct the covariance matrix for the signal, as in Eq. 1 17 1 . 


2.5 Combining model terms 

In order to simplify notation from this point forwards, for each 
pulsar / we combine the matrices G^, and Ff^ into a single, 
Njj X (ntj + 2«sn + 2«dm) matrix, where Njj is the number of TOAs 
in pulsar /, m/ is the number of timing model parameters, and the 
factor 2 in front of both nsN and «dm accounts for the sine and 


cosine terms included for each model frequency. We denote this 
combined matrix T/, such that: 

T, = (G?,Ff,F™), (25) 

and similarily we append the vectors Se,i, asN,/, and aoM,/ to form 
the single vector b/. In this way we can write our complete signal 
model for a single pulsar / as: 

Ti = T,bi. (26) 

We can then construct the block diagonal matrix T such that 
each block is given by the matrix T/ for each pulsar 1, and finally 
append the set of vectors b; for all pulsars to form the complete 
vector of signal coefficients b. In this way the concatenated signal 
model as described thus far for all pulsars, which we denote here 
as T, can be written simply: 

T = Tb. (27) 

2.6 Common noise 

In Tiburzi (2015 PhD Thesis) and Tiburzi et al. (2015, in prep.) it 
was shown that additional sources of noise which are common to 
all pulsars in the PTA can be highly correlated with the quadrupole 
signature of a stochastic GWB. If these sources of noise are present 
in our dataset, we will become less sensitive to a GWB if we do not 
include them in our model. Therefore, in order to ensure that our 
analysis remains robust to the presence of such signals, we will in¬ 
clude in our model the 3 most likely sources of additional common 
noise: 

1: A common, uncorrelated noise term. This allows us to account 
for the possibility that all the millisecond pulsars in our dataset 
suffer from a similar, potentially steep, red noise process, as 
discussed in |Shannon & Cordes|P010l >. 

2: A clock error. [Hobbs et al.| \2Q\2) showed that a PTA is 
sensitive to errors in the time standard used to measure the arrival 
times of pulses. Errors in this time standard would result in a 
monopole signal being present in all pulsars in the dataset. 

3: An error in the Solar System ephemeris. [Champion et al.j 
( [2010^ demonstrated that any error in the planet masses, or any 
unmodelled Solar System bodies will result in an error in our 
determining the barycentric time of arrival of the pulses. This leads 
to a dipole correlation being induced in the timing residuals. 


We note that there are other possible sources of common cor¬ 
related noise in a PTA dataset beyond the three listed above. In 
Section [T7] we will describe models that allow us to fit for a corre¬ 
lated signal, where the form of the correlation is unknown, and is 
described by free parameters in our analysis. In principle one could 
then simultaneously fit for both a GWB, and this additional more 
general signal. While this would significantly decrease our sensitiv¬ 
ity to the GWB it would ensure that our analysis remained robust 
to the existence of unknown correlated signals in the data. More 
optimally, one could perform an evidence comparison between a 
model that includes a GWB, and a model that includes a signal 
with an arbitrary correlation between pulsars in the PTA, in order 
to test which model the data supports. 

A common, uncorrelated noise term can be trivially included 
by adding the model power spectrum to the diagonal of the ele- 
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merits of the matrix that correspond to the intrinsic red noise, 
such that we have: 


H/SN _ „SNc r , „UCr c 

= Vi,i OijOu + iPi 6ij6,j, 


(28) 


where the set of coefficients represent the theoretical power 
spectrum of the common uncorrelated signal, which is the same for 
all pulsars in the array. 

In order to include a clock error within the framework de¬ 
scribed thus far, we append to our matrix T an additional set of 
matrices - one for each pulsar in the array - each of them identical 
to the matrix F™, given by Eq. 1 14 1 , for the corresponding pulsar I. 
Each of these matrices is multiplied by the same set of signal coef¬ 
ficients Ecik, which are appended to the vector of coefficients b, rep¬ 
resenting a single signal being fit to all pulsars simultaneously. We 
use the same number of frequencies in the model for the clock error 
as for the intrinsic spin-noise, and assume a 2-parameter power law 
model for the power spectrum, which we denote as in Eq. 1 16 1 . 
From this we constract the covariance matrix which we define: 


= {afaf) = (29) 

the elements of which can be appended to the total covariance ma¬ 
trix for the signal coefficients *P. We stress that modelling the clock 
signal in this way ensures that we correctly account for both the 
uneven time spans, and unequal weighting of the individual pul¬ 
sars. Additionally, because we fit for the timing model simultane¬ 
ously with the clock signal, the uncertainty in the low frequency 
variations of the signal are factored into the analysis appropriately. 
We show this in a simple simulation in which we use the time 
sampling from our dataset described in Section and include a 
clock error consistent with 10 times the difference between the TAI 
and BIPM2013 time standards, and white noise consistent with the 
TOA uncertainties in our dataset. In Fig. 13 we show the clock 
signal used in our simulation after the maximum likelihood timing 
model has been subtracted from the joint analysis (black line), and 
the time averaged maximum likelihood recovered clock signal with 
Itr uncertainties (red points with error bars). The uncertainties in 
the clock error vary by a factor ~ 9 across the dataset, as differ¬ 
ent pulsars contribute different amounts to the constraints. We find 
the recovered signal is consistent with the injected signal across the 
whole dataspan. 

Finally, in order to model an error in the Solar System 
ephemeris, we can define an error signal e, which will be observed 
in any pulsar I as the dot product between this error vector, and the 
position vector of the pulsar k;, such that the induced residual as a 
function of time, tJ’’’' will be given by: 

t;'’" = e • k,. (30) 


We can incorporate this effect into our analysis by defining a set of 
basis vectors separately for each of the 3 components of e, similarly 
to Eq^] For example, the component in the direction for pulsar 
1 will have basis vectors: 


pcph.x _ cSNt 

r, -r, 

such that the signal induced in the pulsar will be given by: 

eph,x _ ^eph.x 

I ~ * ! %eph,x) ? 


(31) 

(32) 


with a(cph,x) the set of signal coefficients to be fit for. This model 
term is incorporated into our analysis in exactly the same way as 
for the clock error, with the basis vectors for the 3 components 
appended to the total matrix T, the 3 sets of signal coefficients ap¬ 
pended to the vector b, and the diagonal covariance matrix 
constructed from the power law model appended to the matrix ‘P. 



Figure 2. Simulated clock error used used in our analysis (black line) after 
subtracting the maxmium likelihood timing models from the joint analysis, 
and the time averaged maximum likelihood clock signal with Icr uncertain¬ 
ties (red points with error bars). We find the recovered signal is consistent 
with the injected signal across the whole dataspan. 


While this parametrisation does not constitute a physical 
model of the Solar System dynamics, it allows us to incorpo¬ 
rate our uncertainty regarding possible errors in the Solar System 
ephemeris, such as errors in the mass measurements of a number of 
planets or the effects of unknown Solar System bodies. Given the 
dominant source of error in the Solar System ephemeris is likely to 
come from errors in the masses of planets such as Saturn, it could be 
advantageous to include these parameters explicitly in our model. 
In our analysis presented in Section|^we opt for the more conser¬ 
vative approach, and use the general model described here to model 
such errors. 

Once again we include the same number of frequencies in 
the model as for the spin-noise model, and parameterise the power 
spectrum for each of the 3 components, (x,y,z), of the error vector 
e with a separate 2 parameter power law, as in Eq. l|16^. 


2.7 Gravitational-wave background 

When dealing with a signal from a stochastic GWB, it is advanta¬ 
geous to include the cross correlated signal between the pulsars on 
the sky. We do this by using the overlap reduction function - a di¬ 
mensionless function which quantifies the response of the pulsars to 
the stochastic GWB. For isotropic stochastic GWBs, when the pul¬ 
sars are separated from the Earth and from each other by many GW 
wavelengths (i.e., in the short-wavelength approximation, cf|Min-| 
Igarelli & Sidery|2014 1 , this is also known as the Hellings-Downs 
curve jHellings & Downs|1983| l: 


^ [l + +4(1 - cos ^/j) In |sin ^)](1 +Su). (33) 

Here (u is the angle between the pulsars I and J on the sky and 
F(^/j) is the overlap reduction function, which represents the ex¬ 
pected correlation between the TOAs given an isotropic stochastic 
GWB, and the 6ij term accounts for the pulsar term for the autocor¬ 
relation. With this addition, our covariance matrix for the Fourier 
coefficients becomes 
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Vi/SN _ „SN r c , „UC c c , i-/ r \ .GWB c 

^ij.i.j - v>i,i SijSu + ifii SijSjj + r((ij)tpj Sij 


(34) 


both the 2-parameter power law model given in Eq. ( |16) , and also 
take a more general approach, where the power at each frequency 
included in the model is a free parameter in the analysis. In this 
case we define simply as: 




(35) 


where we fit for the set of parameters p, and use a prior that is 
uniform in the amplitude p. 

If we do not want to assume the isotropic (Hellings-Downs) 
overlap reduction function as the description of the correlations be¬ 
tween pulsars in our dataset, we can instead fit for its shape. In 
Section]^ we will do this in two ways: firstly fitting directly for the 
correlation coefficient between each pulsar, r(^;j), and secondly us¬ 
ing a set of four Chebyshev polynomials, where we fit for the coef¬ 
ficients Cl ..4 parameterised such that, defining x = (f/j -7r/2)/(;r/2) 
we will have: 


peaked. Thus, the evidence automatically implements Occam’s ra¬ 
zor: a simpler theory with a compact parameter space will have a 
larger evidence than a more complicated one, unless the latter is 
significantly better at explaining the data. 

The question of model selection between two models “TTo and 
'Hi can then be decided by comparing their respective posterior 
probabilities, given the observed data set D, via the posterior odds 
ratio R: 

^ _ P(H I D) _ PjP I 'Mi)P(Hi) _ Z| P(Hi) 

P{Ho\D) P{D\Ho)P(Ho) ZoPiHo)’ ^ 

•where P{Hi)IP(Ho) is the a priori probability ratio for the two 
models, which can often be set to unity but occasionally requires 
further consideration. 

The posterior odds ratio then allows us to obtain the probabil¬ 
ity of one model compared with the other, simply as: 


P = 


R 

l+R' 


(40) 


3.1.2 MultiNest 


r(x) = Cl + C 2 X + C3(2x^ - 1) + C4(4x^ - 3x). (36) 


3 ANALYSIS METHODS 

While the majority of the results presented in Sectionj^have been 
obtained using a Bayesian approach, we also employ a frequentist 
maximum-likelihood estimator of the GWB strain-spectrum am¬ 
plitude as a consistency check. In the following sections we outline 
the key elements of both these approaches to aid further discussion. 


3.1 Bayesian approach 

3.1.1 General remarks 


Bayesian Inference provides a consistent approach to the estimation 
of a set of parameters 0 in a model or hypothesis H given the data, 
D. Bayes’ theorem states that: 


Pr(0 \D,H) = 


Pr(D I @,‘H)Pr(0 | H) 
Pr(D I H) 


(37) 


where Pr(0 | D,H) = Pr(0) is the posterior probability distribution 
of the parameters, Pr(£) | 0, H) = L(0) is the likelihood, 

Pr(0 I H) = 7r(0) is the prior probability distribution, and 
Pr(£) I H) = Z is the Bayesian Evidence. 

In parameter estimation, the normalizing evidence factor is 
usually ignored, since it is independent of the parameters 0. In¬ 
ferences are therefore obtained by taking samples from the (un¬ 
normalised) posterior using, for example, standard Markov chain 
Monte Carlo (MCMC) sampling methods. 

In contrast to parameter estimation, for model selection the 
evidence takes the central role and is simply the factor required to 
normalise the posterior over 0: 

Z = j L(0);r(0)d"0, (38) 

where n is the dimensionality of the parameter space. 

As the average of the likelihood over the prior, the evidence 
is larger for a model if more of its parameter space is likely and 
smaller for a model where large areas of its parameter space have 
low likelihood values, even if the likelihood function is very highly 


The nested sampling approach jSkilling||2004| > is a Monte-Carlo 
method targeted at the efficient calculation of the evidence, but also 
produces posterior inferences as a by-product. In |Feroz & Hobson| 
( |2008^ and |Eeroz, Hobson & Bridge's! j2009| l this nested sampling 
framework was built upon with the introduction of the MultiNest 
algorithm, which provides an efficient means of sampling from 
posteriors that may contain multiple modes and/or large (curving) 
degeneracies. Since its release MultiNest has been used success¬ 
fully in a wide range of astrophysical problems, from detecting 
the Sunyaev-Zel’dovich effect in galaxy clusters jAMI Consortium! 
|2012[ l, to inferring the properties of a potential stochastic GWB in 
PTA data in a mock data challenge (L13). 

In higher dimensions (> 50), the sampling efficiency of Multi¬ 
Nest begins to decrease significantly. To help alleviate this prob¬ 
lem, MultiNest includes a ‘constant efficiency’ mode, which en¬ 
sures that the sampling efficiency meets some user set target. This, 
however comes at the expense of less accurate evidence values. Re¬ 
cently, the MultiNest algorithm has been updated to include the 
concept of importance nested sampling (INS; |Cameron & Pettitt| 
|20I3|l which provides a solution to this problem. Details can be 
found in |Feroz et al.| ( [20T^ , but the key difference is that, where 
with normal nested sampling the rejected points play no further 
role in the sampling process, INS uses every point sampled to con¬ 
tribute towards the evidence calculation. One outcome of this ap¬ 
proach is that even when running in constant efficiency mode the 
evidence calculated is reliable even in higher (~ 50) dimensional 
problems. In pulsar timing analysis, and especially when determin¬ 
ing the properties of a correlated signal between pulsars, we will 
often have to deal with models that can contain > 40 parameters. 
As such, the ability to run in constant efficiency mode whilst still 
obtaining accurate values for the evidence when these higher di¬ 
mensional problems arise is crucial in order to perform reliable 
model selection. 

All the analyses presented in Section]^ are performed using 
INS, running in constant efficiency mode, with 5000 live points 
and an efficiency of 1%. 


3.1.3 Likelihood function 

Equivalent to the approach described in LI 3, we can write the joint 
probability density of: 
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(i) the linear parameters b, which describe variations in the de- 
terminstic timing model and the signal realisations for the red noise 
and DM variations for each pulsar, and the common noise terms. 

(ii) the stochastic parameters, (a, ft) that describe the intrinsic 
white noise properties for each pulsar, 

(hi) the power-spectrum hyper-parameters that define the spin- 
noise and DM variation power laws, and the spectra of the common 
noise terms such as the stochastic GWB, which we collectively re¬ 
fer to as 0, 

as: 

Pr(b,a,y3,0, I ht) oc Pr(ht|a,/3, b) (41) 

X Pr(b|0) Pr(0)Pr(a,/3)Pr(b). 

In our analysis we simply use priors that are uniform in all 
the model parameters, so we can write the conditional distributions 
that make up Eq. 0 as: 


Pr(ht|a-,/3, b) oc 


1 


Vdet(N) 

X (ht - Tb)], 


exp 


--(<5t-Tb)^N-‘ 


and: 

Pr(b|0) oc 


VdePP 


exp 


—b^T ‘b 

2 


(42) 

(43) 


We can now marginalise over all linear parameters b analyt¬ 
ically in order to find the posterior for the remaining parameters 
alone. 

Defining Z as (T^N^'T -H T-‘), and b as T^N-'ht our 
marginalised posterior for the stochastic parameters a,ft,® alone 
is given by: 


Pr(a,/?, 0|(5t) oc 


det (Z)“ 


Vdet(T) det(N) 

X exp -;^(ht^N^‘ht-b^Z^‘b) 


(44) 


3.2 Frequentist techniques 


As a consistency check of our Bayesian method, we also employ 
a weak-signal regime maximum-likelihood estimator of the GWB 


strain-spectrum amplitude, known as the optimal-statistic i 

Anholm| 

|et al.|2009| [Siemens et al.|2013[|Chamberlin et al.|2014 

1 . It also 


maximises the signal-to-noise ratio (S/N) in this regime, repro¬ 
ducing the results of an optimally-filtered cross-correlation search 
without explicitly introducing a filter function. 


The form of this statistic is 
i:„(5tTp7‘S„P7‘(5tj 

A' = ' % (45) 

2„tr[P7‘S„P7‘S,,] 

where P; = ((5t/(5tJ) is the autocovariance of the post-fit residuals 
in pulsar /, which we can write in terms of the matrices Tj and 
as: 

Pi = Ti'PjTi'^, (46) 

where the matrix *Pi is constructed from maximum-likelihood noise 
estimates obtained in previous single-pulsar analysis. Any GW sig¬ 
nal will have been absorbed into the red-noise estimation during 
this previous analysis. The signal term Sjj is defined such that 
A^Sjj = (St,St]) = S„, where we assume that no signal other than 
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GWs induce cross-correlations between pulsar TOAs. The normal¬ 
isation of is chosen such that (A^) = A^. 

The standard deviation of the statistic in the absence of a 
cross-correlated signal reduces to 


o-Q 


/ 

2tr[P7‘S„P7‘S.i] 

V IJ . 


1/2 


(47) 


which can be used as an approximation to the error on A^ in the 
weak-signal regime. Hence, for a particular signal and noise reali¬ 
sation where we have measured the optimal-statistic, the S/N of the 
power in the cross-correlated signal is given by 


P ^ ZuSt]Vf%jV-/dtj 

(2/ytr[P7‘S„P7‘S,i])'^'’ 


(48) 


with an expectation over all realisations of 

^l/2 

’y/l 


(P>=A2 


2tr[P7‘S„P7'S, 


\ IJ 


(49) 


This S/N effectively measures how likely it is (in terms of num¬ 
ber of standard deviations from zero) that we have found a cross- 
correlated signal in our data rather than an uncorrelated signal. The 
properties of the signal cross-term S;j are determined by a fixed 
input spectral shape, which in this case is a power-law with slope 
7 = 13/3, matching the predicted spectral properties of the strain- 
spectrum resulting from a population of circular GW-driven SMB- 
HBs. 

To compute upper-limits with the optimal-statistic, we follow 
the procedure outlined in |Anholm et al.|j2009^ , where the distribu¬ 
tion of A^ is assumed to be a Gaussian with mean A^ and variance 
<Tq. The latter assumption is clearly only appropriate in the weak- 
signal regime, but serves as a useful approximation. We want to 
find AJi such that, in some predetermined fraction of hypothetical 
experiments (C), the value of the optimal-statistic would exceed the 
actual measured value. Hence we can claim that A^ < A^j to con¬ 
fidence C, otherwise we would have seen it exceed the measured 
value a fraction C of the time. The solution is given by 

A] = A^ + V2o-oerfc-‘ [2(1 - C)]. (50) 


Chamberlin et al. ( |2014^ that the 


Demorest et al.| (|2013|l is identical to the 


It was shown in 
correlation statistic of 
aforementioned optimal-statistic, and in fact allows us to achieve 
a measure of the individual cross-power values between pulsars. In 
the high S/N limit one would expect these cross-power values to 
map out the Hellings and Downs curve when plotted as a function 
of pulsar angular separations. The cross-power values and their as¬ 
sociated errors are given by 


Xu 


6t]Fj%jP-/6tj 

tr[P7‘S„P7‘S,i]’ 


(51) 


0-0,/y = (tr[P7‘SiyP7'Syi])^‘^" 


(52) 


where A ^TuSu = S,j = A%j, and T/j are the Hellings and Downs 
cross-correlation values. 


4 THE DATASET 

Our limits for an isotropic stochastic background are obtained using 
a subset of the full 2015 EPTA data release described in Desvignes 
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Figure 3. Top: timing residuals as a function of Modified Julian Date (MJD) for the 6 pulsars included in the stochastic GWB analysis presented in this 
work, after the maximum likelihood DM variations signal realisation has been subtracted. From top to bottom these are PSRs: J0613-0200, J1012+5307, 
J1600-3053, J1713+0747, J1744-1134, and J1909-3744. While the overall timing baseline for this dataset is ~ 18 yr, only four of the 6 pulsars have data 
that extends across the majority of this timespan, and in particular, PSR J1909-3744 contributes only to the latter half of the dataset, significantly reducing our 
overall sensitivity to signals at the lowest frequencies supported by the dataset. Bottom: Frequency coverage for the 6 pulsars included in the stochastic GWB 
analysis presented in this work. The order of the pulsars is as in the top plot. Colours indicate observing frequencies < lOOOMHz (red crosses), between 1000 
and 2000 MHz (green circles) and > 2000 MHz (blue squares). In addition to fewer pulsars extending across the full dataset, there is also less multi-frequency 
coverage in the early data. This further decreases our sensitivity to a stochastic GWB at the lowest sampled frequencies as the signal becomes highly covariant 
with the DM variations for the individual pulsars in the first half of the dataset. 
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Table 1. Details of the 6 pulsars used for the isotropic stochastic background analysis. Numbers in parentheses correspond to the maximum likelihood values 
from the 5-dimensional analysis described in Section]^ 


Pulsar 

J0613-0200 

J1012±5307 

J1600-3053 

J1713±0747 

J1744-1134 

J1909-3744 

Dataspaii (yr) 

16.05 

16.83 

7.66 

17.66 

17.25 

9.38 

Nsys^ 

14 

15 

4 

14 

9 

3 

o-tps)” 

1.691 

1.610 

0.563 

0.679 

0.801 

0.131 

Logic Asn 

-13.58 ±0.40 (-13.41) 

-13.05 ±0.09 (-13.04) 

-13.71 ±0,54 (-13.42) 

-14.31 ±0.46 (-14.20) 

-13.63 ±0.27 (-13.60) 

-14.22 ±0.42 (-13.98) 

rsN 

2.50 ± 0.99 (2.09) 

1.56 ±0.37 (1.56) 

1.91 ± 1.05 (1.38) 

3.50 ± 1.16(3.51) 

2.21 ±0.82(2.16) 

2.23 ±0.89 (2.17) 

Logic A.DM 

-11.61 ±0.12 (-11.57) 

-12.25 ±0.47 (-11.92) 

-11.75 ±0.39 (-11.67) 

-11.97 ±0.14 (-11.90) 

-12.19±0.38(-11.93) 

-12.76 ±0.53 (-12.51) 

TDM 

1.36 ±0.48 (1.11) 

1.26 ±0.97 (0.27) 

1.64 ±0.80 (1.46) 

2.03 ±0.55 (1.82) 

1.41 ± 1.09 (0.36) 

2.23 ± 1.07 (2.16) 

Global EFAC 

1.01 ±0.02(1.01) 

0.98 ± 0.02 (0.98) 

1.03 ±0.04(1.03) 

1.00 ±0.02 (1.00) 

1.01 ±0.03(1.00) 

1.02 ±0.04 (1.01) 

95% upper limit ^ 

9.7 X 

8.3 X 

2.1 X 10-“* 

4.4 X 10-'^ 

7.0X 10-'5 

5.2 X 


Number of unique observing ‘systems' that make up the dataset for each pulsar. 
Weighted rms for the DM subtracted residuals for each pulsar (D15). 

^ Upper limit obtained from the 5-dimensional analysis described in Sectionj^ 


et al. (in prep.) (henceforth D15). In particular we use a set of 6 pul¬ 
sars, listed in Table[T] that contribute 90% of the total S/N for this 
dataset (see Babak et al. (in prep.) for details). We use this subset of 
the full 42 pulsar dataset in order to minimise the dimensionality of 
the problem, and thus enable accurate evidence calculations using 
MultiNest. The pulsar that contributed next in terms of sensitiv¬ 
ity, PSR J1640+2224, contributes at only the ~ 2% level, so even 
were we to add a small number of additional pulsars, the overall 
gain in sensitivity would be minimal. The DM-subtracted residu¬ 
als, as well as the frequency coverage as a function of time for 
these pulsars are shown in Fig.|^(left and right panel respectively). 
For each of these pulsars a full timing analysis has been performed 
using the TempoNest plugin for the Tempo2 pulsar timing package, 
which simultaneously includes the white noise modifiers EFAC and 
EQUAD for each observing system, as well as intrinsic red noise, 
and frequency dependent DM variations. In Fig. 0 we show the 
mean parameter estimates and Icr uncertainties for the EFAC pa¬ 
rameters obtained for each system from this initial analysis. We 
find all EFACs are consistent with values equal to or greater than 
1 within their uncertainties, with the exception of the Westerbork 
1380 MFlz data in PSR J1713+0747. This could be the result of 
systematic effects that occur in the template forming phase, and is 
the subject of ongoing work. As these systems do not contribute a 
large fraction of the total weight in the dataset, however, it will not 
have a significant impact on the subsequent analysis. Further anal¬ 
ysis of the white noise parameters will be presented in Caballero 
et al. (in prep.). In the joint analysis presented in this work we use 
the linear approximation to the timing model. As such, timing so¬ 
lutions obtained from the initial TempoNest analysis were checked 
for convergence, and the linear regime was found to be suitable in 
all cases. 

The EPTA dataset contains observations from four of the 
largest radio telescopes in Europe: the Effelsberg Radio Telescope 
in Germany, the Lovell Radio Telescope at the Jodrell Bank Obser¬ 
vatory in the UK, the Nancay Radio Telescope in Erance, and the 
Westerbork Synthesis Radio Telescope in The Netherlands. Each 
of these telescopes operates at multiple observing frequencies, and 
so the number of unique ‘systems’ present for any one pulsar can 
be as large as 15. For the 6 pulsars in this dataset we list the number 
of observing systems present in each in Table [T] which in combi¬ 
nation results in 118 white noise parameters. When accounting for 
the 4 spin-noise and DM variation power law parameters for each 
pulsar, we have a total of 142 intrinsic noise parameters before the 
addition of any correlated model components. In an effort to de¬ 
crease the dimensionality we therefore use the mean estimates for 
the EFAC and EQUAD parameters from the individual timing anal- 



Figure 4. EFAC values obtained for all systems from the initial analysis per¬ 
formed for the 6 pulsars used in our analysis. All EFACs are consistent with 
values equal to or greater than 1 within uncertainties, with the exception of 
the Westerbork 1380 MHz data in PSR J17I3+0747 which have values con¬ 
sistent with ~ 0.5. This could be the result of systematic effects that occur 
in the template forming phase, and is the subject of ongoing work. 


ysis presented in D15, and fit a single global EFAC for each pulsar, 
reducing the number of intrinsic noise parameters to 30. 

In order to check the validity of this simplification we per¬ 
formed a 5-dimensional analysis for each of the 6 pulsars in this 
dataset, fitting for power law intrinsic red noise and DM variations, 
in addition to a global EFAC parameter after adjusting the error bars 
using the mean values from D15. The parameter estimates obtained 
are given in Table [T] and the 1-dimensional marginalised posteriors 
for J1909-3744, J1713+0747, and J1744-1134 from this analysis 
are shown in Fig. In each case we show the red noise and DM 
variation power law parameters for the full noise analysis (red) and 
5-dimensional analysis (blue). We also show the global EFAC pa¬ 
rameter from the 5-dimensional analysis in each case. We find the 
posteriors are consistent between the two sets of analysis. 

Table [T]also lists the 95% upper limit on a red noise pro¬ 
cess with a spectral index of 13/3 at a reference frequency of lyr^’ 
for each of the 6 pulsars used in our analysis. This limit was ob¬ 
tained when simultaneously fitting for the 5 intrinsic noise param¬ 
eters for each pulsar in addition to the steep spectrum noise term. 
The two pulsars with the most constraining upper limit are PSRs 
J1909-3744, and J1713+0747, consistent with the results obtained 
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Figure 5. Comparison of the 1 -dimensional marginalised posterior probability distributions for PSRs (left to right) J1909-3744, J1713+0747, and J1744-1134. 
The y-axis in all plots represents probability. In each case we show the spin-noise and DM variation power law parameters for the full noise analysis (red solid 
lines) and 5-dimensional analysis where the TOA error bars have been pre scaled by the mean value of the EFAC/EQUAD parameters for each pulsar backend 
(blue dashed lines). In both cases pai'ameter estimates have been obtained using a uniform prior on the amplitude of the spin-noise and DM variations power 
law models. We also show the global EFAC pai'ameters from the 5-dimensional analysis in each case. We find the posteriors are consistent between the two 
sets of analysis. 
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Table 3.95% upper limits on the amplitude of an isotropic stochastic GWB 
obtained for different models at a reference frequency of lyr“^. 


Model 95% upper limit 

(xlO^‘^) 

Bayesian Analysis 


Fixed Noise - Fixed Spectral Index 1.7 

Varying Noise - Fixed Spectral Index 3.0 

Additional Common Signals - Fixed Spectral Index 3.0 


Fixed Noise - Varying Spectral Index 8.0 

Varying Noise - Varying Spectral Index 13 

Additional Common Signals - Varying Spectral Index 13 


Frequentist Analysis 

Fixed Noise - Fixed Spectral Index 2.1 


Simulations - Varying Spectral Index 

White Noise Only 4.3 

White and intrinsic spin-noise 7.2 

White and intrinsic spin-noise and DM variations 12 


Table 4. 95% upper limits obtained for common noise terms at a reference 
frequency of lyr“*. 

Model 

95% upper limit 
(xl0-‘^) 

Additional Common Signals - Varying 

Spectral Index 

Acn 

13 

Aclk 

11 

■^eph (x) 

65 

-^eph (y) 

14 

■^eph (z) 

25 


in Babak et al. (in prep.), with values of 5 x 10 and « 4x 10 
respectively. 


5 RESULTS 

5.1 Limits on an Isotropic Stochastic GWB 

5.1.1 Bayesian approach 

In Table 1^ we list the complete set of free parameters that we in¬ 
clude in the different models used in the analysis presented in this 
section, along with the prior ranges used for those parameters. 

When parameterising the stochastic GWB using the power law 
model in Eq. l|^, we run two parallel sets of analyses: in the first set 
we fix 7 = 13/3, consistent with a stochastic GWB dominated by 
SMBHBs; in the second set we allow 7 to vary freely within a prior 
range of [0,7], In both cases we consider three different models: 

(i) with the intrinsic timing noise for each pulsar fixed at the 
maximum likelihood values given in Table[T] 

(ii) with the intrinsic timing noise for each pulsar allowed to 
vary; 

(iii) as in (ii), but including additional common uncorrelated red 



Figure 7. One and two-dimensional marginalised posterior parameter es¬ 
timates for the amplitude and spectral index of a correlated GWB in the 
6 pulsar dataset presented in this paper when varying the intrinsic noise 
parameters for each pulsar. The amplitude and spectral index are highly 
correlated, resulting in a significantly higher upper limit when allowing the 
spectral index to vary, as opposed to fixing it at a value of 13/3. 


noise, a clock error, and errors in the Solar System ephemeris as 
discussed in SectionE^l 


The 95% upper limits for the amplitude of an isotropic 
stochastic GWB in the six different models are listed in Table |3 
In Table|^we then list the 95% upper limits for the additional com¬ 
mon noise terms that were included in model (iii), when allowing 
the spectral indicies to vary. All upper limits in this section are re¬ 
ported at a reference frequency of lyr“*. 

The one-dimensional marginalised posteriors for the ampli¬ 
tude of the GWB for each of these models are shown in Fig. 

We find that in both the fixed, and varying spectral index model 
for the GWB, limits placed under the assumption of fixed intrinsic 
timing noise are erroneously more stringent than when the noise is 
allowed to vary, by a factor ~ 1.8 and ~ 1.6 respectively. This is 
a direct result of using values for the intrinsic noise that have been 
obtained from analysis of the individual pulsars, in which the red 
spin-noise signal, and any potential GWB signal will be completely 
covariant. The natural consequence is that fixing the properties of 
the intrinsic noise to those obtained from the single pulsar analysis 
will always push the upper limit for the GWB lower in a subsequent 
joint analysis. 

Both of the most recent isotropic GWB limits from pulsar tim¬ 
ing have been set using frequentist techniques, either by performing 
a fixed noise analysis jPemorest et al.|2013| >, or using simulations 
I Shannon et al. 2013 1 , obtaining 95% upper limits of 7 x 10“*^ and 
2.4 X 10“'^ respectively. In both cases, therefore, the analysis per¬ 
formed was fundamentally different to the Bayesian approach pre¬ 
sented in this work. As such it is difficult to compare our results 
directly, or to ascertain the effect of fixing the intrinsic noise pa¬ 
rameters on limits obtained using those methods. 

The most recent limit placed when allowing the intrinsic noise 
parameters of the pulsars to vary is given by |van Haasteren et al.| 
( 201l| , in which a 95% upper limit of A = 6 x 10“'^ was ob¬ 
tained at a spectral index of 13/3. Our model (ii) is most com- 
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Table 2. Free parameters and prior ranges used in the Bayesian analysis. 


Pai'ameter 

Description 

Prior range 


White noise 


Global EFAC 

uniform in [0.5 , 1.5] 

1 parameter per pulsar (total 6) 

Spin-noise 


Spin-noise power law amplitude 

uniform m[10-20 , iQ-'®] 

1 parameter per pulsar (total 6) 

ysN 

Spin-noise power law spectral index 

uniform in [0,7] 

1 parameter per pulsar (total 6) 

DM variations 


DM variations power law amplitude 

imiform in [10--“ , 

1 parameter per pulsar (total 6) 

yoM 

DM variations power law spectral index 

uniform in [0,7] 

1 parameter per pulsar (total 6) 

Common noise 


Uncorrelated common noise power law amplitude 

uniform in [10--“ , 

1 parameter for the array 

ycN 

Uncorrelated common noise power law spectral index 

uniform in [0,7] 

1 parameter for the array 

'^clk 

Clock error power law amplitude 

uniform in [10-^° , 

1 parameter for the array 

yelk 

Clock error power law spectral index 

uniform in [0,7] 

1 parameter for the array 

'^eph 

Solar System ephemeris error power law amplitude 

uniform in [10-^° , IQ-'"] 

3 parameters for the array (x, y, z) 

y^ph 

Solai- System ephemeris error power law spectral index 

uniform in [0,7] 

3 pai'ameters for the array (x, y, z) 

Stochastic GWB 

A 

GWB power law amplitude 

uniform in [10--“ , 

1 parameter for the array 

y 

GWB power law spectral index 

uniform in [0,7] 

1 parameter for the array 

Pi 

GWB power spectrum coefficient at frequency i/T 

uniform in [10-^° , 10“] 

1 parameter for the array per frequency in 
unparameterised GWB power spectrum model (total 20) 

Stochastic background angular correlation function 

C1...4 

Chebyshev polynomial coefficient 

uniform in [-1,1] 

see Eq. j36j 

Tij 

Correlation coefficient between pulsars (I,J) 

uniform in [-1,1] 

1 parameter for the array per unique pulsar pair (total 15) 




5e-15 1e-14 1.5e-14 2e-14 

GWB Power Law Amplitude 


Figure 6. One-dimensional marginalised posterior parameter estimates for the amplitude of a correlated GWB in the 6 pulsar dataset presented in this paper 
when: (left) Fixing the spectral index of the power law to a value of 13/3, consistent with a background dominated hy a population of SMBHB, and (right) when 
marginalising over the spectral index given a prior range of [0,7], In each case we show the posterior given: (red solid line) Fixed intrinsic noise parameters for 
each pulsar, where the values of the parameters are given by the maximum likelihood estimates listed in Table[T] (black dotted line) varying noise parameters 
for each pulsar, and (magenta dashed line) varying noise parameters for each pulsar, and additionally including a common uncorrelated red noise process, 
clock en'ors, and en'ors in the Solar System ephemeris in the model. Vertical lines in each case represent the 95% upper limits for each model. 


parable to this analysis, in which we obtain a 95% upper limit of 
A = 3 X 10“'^, an improvement of a factor of two. This translates 
into a limit on f2g„(/)/i^ = 1.1 x 10“^ at 2.8 nHz. We confirm this 
result by analysing the 2015 EPTA dataset with model (ii) using an 
independent cod^ which makes use of the PAL, parallel-tempered 
adaptive MCMC sample^which explores the parameter space in a 
fundamentally different way to MultiNest, and obtain a consistent 
95% upper limit. 

Finally for model (iii) when including additional common or 

^ https://github.com/stevertaylor/NX01 
^ https://github.com/jellisl8/PAL2 


correlated terms in the analysis we find the extra parameters have 
a negligible impact on our sensitivity, with consistent upper limits 
obtained in both the fixed and varying spectral index models. 

We find the upper limits for the uncorrelated common red 
noise model to be consistent with those obtained for the GWB, 
however we find the upper limit for the clock error signal to 
be slightly lower, with A^k < 1.1 x 10“*'* compared to Acn < 
1.3 X 10“''*. This is to be expected however, as the clock is han¬ 
dled coherently across all pulsars, whereas the GWB and common 
uncorrelated red noise signals are handled incoherently, as such we 
have greater sensitivity when searching for the clock signal and ob¬ 
tain a correspondingly lower limit for the amplitude. 
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Figure 8. One and two-dimensional marginalised posterior parameter esti¬ 
mates for the amplitude and spectral index of a correlated GWB in 3 sim¬ 
ulated datasets, including (Top) White noise only, (Middle) White and in¬ 
trinsic spin-noise only, and (Bottom) White noise, intrinsic spin-noise, and 
dispersion measure variations. In each case we use the TOAs from the 6 
pulsars used in the GWB analysis presented in this work, and use the maxi¬ 
mum likelihood noise parameters from Table^when constructing the sim¬ 
ulations. The upper limits obtained in each case are 4.3 X 10“’^, 7.2x 10“*^, 
and 1.2 X 10“*'*. We find both the upper limit, and the form of the posterior 
to be consistent between the third simulation and the real dataset. In both 
cases we are simply recovering our uniform prior on the amplitude of the 
GWB signal at small amplitudes, before the data begins to place constraints 
on the upper limit at large amplitudes. 
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The limits on the different errors originating in the Solar Sys¬ 
tem emphemeris can he understood given the components of the 
unit vector from the SSB toward the two pulsars that contribute 
most to our analysis, PSRs J1713-I-0747 and J1909-3744, which 
are given by (-0.20,-0.97,4-0.14) and (4-0.24,-0.75,-0.61) re¬ 
spectively. Both PSRs J1909-3744 and J17134-0747 contribute 
very little to the constraints on the ephemeris in the x direction, and 
so here we see the greatest degradation in the limit on the ampli¬ 
tude, while in the y direction PSR J17134-0747 contributes almost 
fully and so the limit we obtain is only slightly worse than that 
obtained for the GWB and uncorrelated common noise terms. 

We consider model (iii) to be the most robust analysis pre¬ 
sented in this paper, and so conclude that the 95% upper limit pro¬ 
vided by our dataset on a power law GWB is A < 3.0 x 10“*^ at 
7 = 13/3 , and A < 1.3 x 10“*'* when marginalising over spectral 
index. 

That the upper limit is considerably higher in the varying spec¬ 
tral index model can be understood from the two-dimensional pos¬ 
terior distribution for A and y in Fig.|^ Here we see the clear cor¬ 
relation between the two quantities; as we will see below, our PTA 
is most sensitive at frequencies « 1 yr“*, meaning that for a single 
power-law spectrum, the flatter the spectral index, the less stringent 
the limit on A. 

As a consistency check on this result we perform a set of three 
simulations using the sampled time stamps present in the actual 6 
pulsar dataset. In the first simulation we include only a white noise 
component with an amplitude determined using the TOA uncer- 
tanties from the real dataset. In the second simulation we then add 
an intrinsic spin-noise component, with amplitudes and spectral in¬ 
dices equal to the maximum likelihood values presented in Table 
[T] Finally in the third simulation we also include DM variations, 
where as with the intrinsic spin-noise we use the maximum likeli¬ 
hood values in Table to set the amplitudes and spectral indices. 
Critically in all the simulations we include no correlated GWB 
term, and so in each case we expect to recover only our uniform 
prior on the amplitude of the GWB included in our model. 

The one and two-dimensional marginalised posterior parame¬ 
ter estimates for the amplitude and spectral index of the GWB from 
each of the three simulations are shown in Fig.l^ We obtain upper 
limits of 4.3 x 10“*^, 7.2 x 10“*^, and 1.2 x in each case re¬ 
spectively. In all the one-dimensional posterior distributions for the 
amplitude of the signal we are simply recovering our prior, such 
that the probability is proportional to the amplitude of the signal 
below some limit set by the data. In the case of the third simula¬ 
tion, which is most similar to the real dataset, both the upper limit, 
and the form of the posterior is consistent with the results presented 
in Fig. 1^ and Table 

In Fig. 1^ we show the 95% upper limits from a power spec¬ 
trum analysis that does not assume a power law model, but allows 
the power at each frequency included in the model to vary sep¬ 
arately. We perform this analysis separately both for a correlated 
GWB (red points), and uncorrelated common red noise process 
(blue points) while varying the intrinsic noise parameters for each 
pulsar, however we do not include any additional common terms. 
The upper limit obtained from the equivalent power law analysis 
of the GWB at a spectral index of 13/3 of 3 x 10“*^ is overplot¬ 
ted as a straight line. Frequencies were included from l/T up to 
20/r, with T = 17.66 yr, beyond which we do not expect the data 
to provide significant constraints on a steep red noise process. We 
find that our limit at a spectral index of 13/3 is most heavily con¬ 
strained by the 3/7’ term, corresponding to / »: 5 x 10“®Hz. This is 
likely a combination of the lack of multifrequency data in the early 
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Figure 12. (left) Frequentist upper-limits on the strain amplitude of the stochastic GWB obtained via the optimal-statistic for our 6 pulsar dataset. Red lines 
indicate 90% and 95% upper-limits at the fiducial slope of the strain-spectrum of -2/3, which coiresponds to a slope of the residual PSD of -13/3. (right) 
The individual cross-power values are shown for our 6 pulsar dataset. All values are consistent with zero cross-correlation. 



Figure 9. (Top) 95% upper limits from an unparameterised power spectrum 
analysis for a correlated GWB (red points), and uncorrelated common red 
noise process (blue points) for the 6 pulsar dataset described in Section]^ 
obtained while varying the intrinsic noise parameters for each pulsar. The 
upper limit obtained from a power law analysis of the GWB at a spectral 
index of 13/3 of 3 X 10“'^ is overplotted as a straight line. Frequencies 
were included from l/T up to 20/r, with T = 17.66 yr. Beyond these 
frequencies the data provides increasingly poorer constraints on a steep red 
noise process, and so we do not consider higher frequency terms in our 
model. Both the correlated and uncorrelated power spectrum are completely 
consistent with one another at all frequencies. The difference in the log 
Evidence between the correlated and uncomelated models was 0.2 ± 0.3, 
indicating no support for the presence of a correlated signal in the dataset. 


data for PSR J1713-1-0747 as shown in Fig.[^ which significantly 
impacts our ability to disentangle DM variations from frequency- 
independent red noise at the lowest frequencies, and that our PSR 
J1909-3744 dataset is only ~ 9 yr in length. Despite these limi¬ 
tations, the long timing baseline of the EPTA dataset used in this 
work is still critical for placing limits on the lowest frequencies in 
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Figure 10. One-dimensional marginalised posterior parameter estimates for 
the amplitude of a correlated GWB (red solid line) and an uncorrelated com¬ 
mon red noise model (blue dashed line) in the 6 pulsar' dataset presented in 
this paper when marginalising over the spectral index given a prior range 
of [0,7]. Posteriors were obtained when varying the intrinsic noise param¬ 
eters, and including only either the GWB, or common uncorrelated terms. 
We find the upper limits to be completely consistent with one another, and 
obtain a change in the log evidence of -1.0 ± 0.5 for the GWB model over 
the uncorrelated common red noise model, suggesting no strong support for 
either. 


our analysis. Both the correlated and uncorrelated power spectrum 
are completely consistent with one another at all frequencies. The 
difference in the log Evidence between the correlated and uncorre¬ 
lated models was 0.2 ± 0.3, indicating no support for the presence 
of a correlated signal in the dataset. 

In Fig. [T^ we further assess the impact of including, or not, 
the Hellings-Downs correlation on the upper limit obtained in our 
power law model (ii) for the varying spectral index case. We show 
the one-dimensional marginalised posterior for both the amplitude 
of the GWB, which includes the correlation between pulsars (red 
solid line), and for the amplitude of an uncorrelated common red 
noise power law process (blue dashed line). We find the upper limits 
to be completely consistent with one another, and obtain a change 
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Figure 11. The recovered conelation between pulsar s as a function of an¬ 
gular separation on the sky for a power law noise process. The red and blue 
lines represent the 68% and 95% confidence intervals for the correlation 
function when modelled by the lowest 4 Chebyshev polynomials, while 
the individual points are the mean correlation coefficient with Icr uncer¬ 
tainty for each pulsar' pair when fitting without assuming a smooth model. 
The Hellings-Downs correlation is represented by the dotted line. For both 
models the correlation of a common power law model between pulsars is 
consistent with effectively all possible values (the range [-1,1]). 


in the log evidence of -1.0±0.5 for the GWB model over the uncor¬ 
related common red noise model, indicating no strong support for 
either model in the data. We confirm this result by obtaining con¬ 
straints on the correlation between pulsars as a function of angular 
separation for a common power law model, where the amplitude 
and spectral index of the power law are free to vary, and we fit si¬ 
multaneously for the intrinsic noise parameters for the individual 
pulsars. We fit for the correlation using the two methods described 
in Section [TTt] using either the four lowest order Chebyshev poly¬ 
nomials, or fitting for the correlation coefficient directly. In Fig. 
[m the red and blue lines represent the 68% and 95% confidence 
intervals for the correlation function when modelled by the low¬ 
est 4 Chebyshev polynomials, while the individual points are the 
mean correlation coefficient with Icr uncertainty for each pulsar 
pair when fitting directly. In both cases the correlation is consistent 
with effectively all possible values (the range [-1,1]). 


5.1.2 Frequentist approach 

Applying the optimal-statistic introduced in Section [7!2| to our re¬ 
duced 6 pulsar dataset, and testing for a strain-spectrum slope of 
-2/3, gives = (-2.86 ± 4.29) x 10“^°, with an associated 
S/N = -0.67. This is clearly a non-detection, however the 95% 
upper-limit on A is 2.05 x 10“'^, which is more constraining than 
the best published limit of |Shannon et al.|P013[ l, and is consistent 
with the fixed noise Bayesian limit of 1.7 x 10“'^. The optimal- 
statistic limits as a function of slope are shown in the left panel of 
Fig.[^ with limits at the fiducial slope value marked in red. 

The computed cross-power values for our 6 pulsar dataset are 
shown in the right panel of Fig.[^ and are all consistent with zero 
correlation, as expected. 
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6 DISCUSSION 


6.1 Implications for SMBHB astrophysics 

As diseased previously, the most promising astrophysical source of 
CWs in the nHz regime relevant to PTA observations is a cosmo¬ 
logical population of adiabatically in-spiralling SMBHBs. 

In this section we will consider the implications of the upper- 
limit obtained in Section on the gravitational-wave signal for 
models of astrophysical populations of SMBHBs. If we assume 
that a SMBHB evolves purely due to gravitational radiation reac¬ 
tion, and that all SMBHBs are in circular orbits - we will return 
to these assumptions at the end of the section - the characteristic 
amplitude, Eq. 0, is given by: 
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A limit on the amplitude A, as given by Eq. ( |53^ therefore places 
constraints on d^n/(dzdM), i.e., the number density of SMBHB 
mergers per unit redshift and unit chirp mass across cosmic his¬ 
tory. Although the dominant contribution to the signal comes from 
relatively massive (M > 10® Mq), low redshift (z < 2) systems 
( [Sesana, Vecchio & Colacino|2008^ , their merger rate is still poorly 
constrained, resulting in a fairly wide range of possible signal 
amplitudes. [Sesana, Vecchio & Colacino| ( 2008|( exploited semi- 
analytical merger trees from [Volonteri, Haardt & Madau| j2003| l 
to estimate a plausible range for the amplitude of a CWB of 
5 X 10“'® < A < 3 X 10“'®. Merger rates extracted from cosmologi¬ 
cal simulations like the Millennium jSpringel et al.|2005|( and Mas¬ 
sive Black ( [Khandai et al.|2014| l simulations, coupled with ditferent 
prescriptions for the SMBH-galaxy relation results in a computable 
range of 4 x 10“'® < A < 2 x 10“'® 

[20091 |Ravi et al.|2012^ . Recently, jS 


I Sesa na, Vecchio & Volonteri 
1 2013b^ constrained the 


expected range of A by building a set of phenomenological models 
based on the observed properties of interacting galaxies. Thousands 
of models fulfilling all relevant observational constraints were as¬ 
sembled by combining different estimates of the galaxy mass func¬ 
tion, pair fractions, estimated merger times and galaxy-SMBH re¬ 
lations. The central 90% of the probability distribution function 
(PDE) in the amplitude A lies in the range 3x10“'® < A < 3x 10“'®, 
as shown in the right panel of Fig.[T3] A similar approach was em¬ 
ployed by [Ravi et al.| ( |2014^ , yielding consistent results. 

In Eig.|13|we compare the expected range in predicted by 
the phenomenological models presented in jSesanaj j2013b| l to the 
95% upper limit obtained in Section HU Shaded areas represent 
the central 68%, 95%, 99.7% and 100% confidence intervals for 
the predicted signal. The red curve is derived by converting the 
95% limit on the unparametrised power spectrum shown in Eig.j^ 
into K as 


h, = (power X 127T^fTy'^, 


(55) 


where T is the total observation time. Eq. |55) can be calculated 
directly from Eq. by noting that the power at each frequency 
is the integral of S (/) over a frequency bin A/ = l/T. Fig. [T^ in¬ 
structively shows how the limit on A that is usually quoted in the 
literature is extrapolated from the actual sensitivity of the PTA. Our 
dataset is most sensitive at / »: 5 x 10“®Hz, where the 95% limit on 
hc 'K~ 1.1 X 10“'"'. This is then extrapolated to / = lyr“' assuming 
a power-law, to get a 95% upper limit of A = 3.0 x 10“'®. 
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The right panel of Fig. shows the probability distribution of A 
inferred from the theoretical models of [Sesan^ j2013b| l together 
with the region excluded at 95% confidence by our analysis. Only 
about 5% of the distribution is excluded, meaning that our limit 
does not place severe restrictions on the cosmic SMBHB popu¬ 
lation. Note that A values obtained from merger trees and cos¬ 
mological simulations quoted above are generally in the range 
4 X 10“*^ < A < 3 X 10“'^, and therefore also consistent with this 


limit. 

We caution that Fig. [T^ shows the expected range of the GW 
signal given circular GW driven SMBFlBs, with negligible con¬ 


pling to the environment. It has been shown (see, e.g. 

Enoki & Na- 

gashima|2007||Kocsis & Sesana|20II |Sesana|20I3a| 

McWilliams, 

Ostriker & Pretorius||20I4 Ravi et al.||20I4bb that both high ec- 


centricities and strong environmental coupling might cause a sig¬ 
nificant suppression of the GW signal at / < lO^^Hz. If this is 
the case, our 1.1 x 10“*"* limit at / « 5 x lO^^Hz cannot be ex¬ 
trapolated to / = lyr^' using a simple power-law model. Even a 
much more stringent limit obtained at such low frequencies might 
not have significant implications for the cosmic SMBHB merger 
rate, because the suppression at low frequencies would invalidate 


the / extrapolation to higher frequencies (see Fig. 2 in Sesana 


|2013a^ . Since the detailed dynamical evolution of SMBHBs cou¬ 
pled with their environment is still poorly understood, current PTA 
limits do not allow us to formulate strong astrophysical statements 
about the cosmological population of SMBHBs. 



Observed GW Frequency [Hz] 


6.2 Limits on the cosmic (super)string tension 

The limits computed in Section can be converted into upper 
limits on the linear energy density of a cosmic (super)string net¬ 
work, fi, or tension in the Nambu-Goto approximation, usually de¬ 
scribed by the dimensionless quantity G^lc^, where G is New¬ 
ton’s constant and c the speed of light. Field theory cosmic strings 
( |Kibble|1976| |Jeannerot, Rocher & Sakellariadou|2003^ , are one¬ 
dimensional topological defects, relics of an early, more symmetric 
state of the Universe, created through the mechanism of sponta¬ 
neous symmetry breaking during the various phase transitions that 
the early Universe underwent. Their formation is a generic property 
of supersymmetric hybrid inflation scenarios ([J eannerot, Rocher &| 
|Sakellariadou|20d3) , whereas the creation of their superstring the¬ 
ory counterparts, usually referred to as cosmic superstrings, are also 
a natural by-product of brane inflation scenarios (e.g. [Saran g i ^ 
|Tye|2002{|Jones, Stoica & Tye|2003) . 

A cosmic string network consists of “infinite” (larger than the 
particle horizon) strings and cosmic string loops. A cosmic string 
network grows along the expansion of the Universe and is expected 
to settle in a scaling regime, were all the fundamental properties 
of the network grow proportionally with cosmic time, something 
achieved with the creation of loops; the main energy loss mecha¬ 
nism of the network. When two cosmic strings intersect, they “in¬ 
tercommute” (exchange partners) with a characteristic probability, 
and form loops. Cosmic string loops, once created, oscillate and 
decay emitting all of their energy in various forms of radiation, 
with the dominant form thought to be GWs ^Vilenkin|1981| >. The 
stochastic GWB created by a cosmic string network is broadband 
(from -10“’® Hz, to higher than 10^ Hz, depending on the size of the 
loops created), a characteristic feature of primordial GW sources, 
and is potentially detectable by any present or future GW detector. 
The cosmic string GW spectrum consists of a flat part at high fre¬ 
quencies, originating from loops decaying in the radiation era, and 
a broad peak at lower frequencies originating from loops decaying 


Figure 13. Comparison between the expected GWB amplitude from a cos¬ 
mological population of SMBHBs and the 95% upper limit obtained with 
our PTA experiment. Shaded areas represent the central 68%, 95%, 99.7% 
and 100% confidence interval of the predicted signal according to|Sesana| 
j2013b) , whereas the red curve is the 95% upper limit presented in this 
paper, obtained by converting the unparametrised power spectrum shown 
in Fig. I^to characteristic amplitude. The black triangle is the extrapolated 
95% upper limit on A at / = lyr“*. The right panel shows the PDF of the 
predicted aX f = lyr^', and the shaded area marks the region excluded at 
95% confidence by our limit (less than 5% of the distribution). 


in the matter era. PTAs are in the privileged position to typically 
probe the GW emission originating from loops decaying either in 
the matter era or in the radiation-to-matter era transition, making 
them excellent instruments to detect the stochastic GWB of a cos¬ 
mic string network. In the case of a non-detection, the upper limits 
on the GWB can be used to constrain the energy scale of cosmic 
strings. 

There are nume rous investigations on the cosmic string GWB 
in the literature (e.g. Damour & Vilenkin|2005 Olmez, Mandic & 


|Siemens|2010[|Depies & Hogan^O" l, depending 


: on various as¬ 


sumptions concerning the GW emission mechanisms and the com¬ 
putation of the loop number density. The cosmic string GWB spec¬ 
tra used in this paper are those computed in |Sanidas, Bat t ye & Stap-| 
|pers|(|2013|l, which are based on an updated version of the mod¬ 
elling presented in |Sanidas, Battye & StappersH2012) l. The main 
difference of these two investigations is the inclusion of the effects 
of massive particle annihilation, which reduces the amplitude of 
the spectrum at higher frequencies. Both investigations, however, 
do not make any assumptions about the values of the fundamen¬ 
tal model parameters used to compute the GW spectrum, grant¬ 
ing characteristic robustness to the results presented here. |Sanidas,| 
[Battye & Stappers| ( |2012^ presented a generic way to compute the 
GWB based on the widely accepted one-scale model for cosmic 
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Q. spectral slope 


Figure 14. The 95% confidence upper limits on the amplitude of an 
isotropic, stochastic GWB as a function of the local spectral index at the 
frequency of 1 yr“*, imposed by models (ii) (dashed black) and (iii) (solid 
red) of Section [TT] The limits are expressed in terms of the strain of the 
GWB, he, and the dimensionless spectral energy density of GWs, Clh^, 
where Hq = 100/tkms^* Mpe^*, as a function of their respective spectral 
indices. 


strings. The basic parameters of this model are the string tension 
Gjilc?', the birth-scale of loops relative to the horizon, and 
the intercommutation probability, p. Based on this modelling, and 
assuming that the cosmic string network maintains a scaling evo¬ 
lution along the expansion of the Universe, one can compute the 
loop number density for all the parameter combinations that create 
a GWB at the frequency regime probed by PTAs. Additionally, no 
assumption is made concerning the dominant GW emission mecha¬ 
nism from the strings (i.e., kinks or cusps), which is modelled using 
two additional parameters, a spectral index, q, and a cut-off, n,, on 
the number of emission harmonics n. 

In this section we present updated constraints on the cosmic 
string tension, for various values of the intercommutation proba¬ 
bility, but independent of all the rest model parameters. The GWB 
upper limits we used are those obtained in the analysis for the mod¬ 
els (ii) and (iii) of Section mi where the spectral index y was a 
free parameter, as these are presented in Fig.[T^ Additionally, we 
used the GW sensitivity curve in Fig.[T^to investigate tension con¬ 
straints across the probed frequency range. In the first case, we used 
the information of both the amplitude and the local spectral index 
of the GWB in order to create the tension exclusion curves in the 
cosmic string model parameter space, whereas in the second case 
we used only the amplitude information. In Fig. [T^ we present 
the cosmic string tension exclusion curves in the parameter space 
Gplc^ — ttes accessed by the PTAs, for field theory strings (p = 1). 
These exclusion curves were constructed from the combination of 
the «, = 1 and «, = 10'* exclusion curves, depending on which 
one provided the highest tension value. The limits are provided by 
the n, = \ networks in the mid-ttes region (-8 :£ logio ®c.s “3), 
whereas in the large (logjQ <: -3) and small (-8 < logm ttes) 
loop regimes, by the n, = 10"* networks. As discussed in Sanidas, 


Battye & Stappers (2012 i, the n, = I and n, = lO'^ cases will al¬ 


ways provide the largest tension values for fixed values of the rest 
cosmic string model parameters. The limits from the n, = 1 net¬ 
works are independent of the GW emission mechanism since the 
power emitted per emission mode is oc n^‘i. For the limits provided 
by the n, = 10** networks, we used a spectral index q = 4/3, cor- 



Figure IS. Exclusion limits for different cosmic string network configura¬ 
tions with the same Clh^ value at a frequency / = I yr“* in the Gpje^ — acs 
parameter space, for p = \. Both the amplitude and spectral slope infor¬ 
mation of the GWB limits were used to construct the limits. The dashed 
black and solid red curves, are based on the results of models (ii) and (iii) 
presented in Fig. [irrespectively. In the mid-ffes region, the tension upper 
limits are provided by the n, = 1 networks, whereas in the regime of large 
and small loops, by the n, = 10“* networks. 


responding to cusp emission, which always provides larger tension 
values than the q = 2 case which corresponds to emission from 
kinks. The cosmic string tension upper limit for model (ii), where 
the intrinsic timing noise of each pulsar is allowed to vary, is 

Gpje^ < 1.2 X 10“^(95% confidence), (56) 

whereas the tension upper limit for model (iii) where the common 
uncorrelated red noise, clock and Solar System ephemeris errors 
are included is 

Gpje^ < 1.3 X 10^^(95% confidence). (57) 


For the particular case of large loop production with a = 0.05, 
as suggested by the most recent Nambu-Goto cosmic string evolu¬ 
tion simulations ( |Blanco-Pillado, Plum & Shlaer|201l||2014) , we 
get an upper limit Gyu/c^ < 3.0 x 10“". The tension constraint that 
was obtained based on that loop number density and the previous 
EPTA GWB limit was Gpje^ < 2.8 x 10“^ (jBlanco-Pillado, Olum 


|& Shlaer|2014[ l. Assuming a -1.7 times improvement on the value 
of the Gg* used in that work, one would expect approximately a 
constraint Gp/c^ < 9.7 x 10“*°. This order of magnitude differ¬ 
ence stems mostly from the normalization imposed on the produced 
number of loops. Whereas in our model we assume that all of the 
energy lost by the cosmic string network in order to attain scal¬ 
ing is channeled in loops of just one size, the loop production in 
( |Blanco-Pillado, Olum & Shlaer|2014] l takes place on a much wider 
range of large loops (wider range of a). The difference in these two 
constraints is expected, since as we have demonstrated in [Sanidas^ 
[Battye & Stapperij {2012\ , where we also assumed a log-normal 
distribution for the loop birth scale, a wider range for the values of 
a has as an effect the lowering of the matter era peak GWB ampli¬ 
tude and its broadening. However, a direct comparison of these two 
results is not straightforward, not only because of the loop num- 
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Table 5. Upper limits on the cosmic superstring tension Gjii/c- for p + \ 
and the two scaling laws proposed in the literature, using the GWB limits 
placed by models (ii) and (iii). 


Model 

Scenaiio ii 

(varying spectral index, 
varying noise) 

Scenario iii 

(varying spectral index, 
additional common noise) 

Scaling law 

k=0.6 

k=l 

k=0.6 

k=l 

p = 10“' 

2.2 X 10“® 

1.1 X 10“" 

2.4 X 10-* 

1.0 X 10-" 

p = 10-2 

7.3 X lO-** 

1.6 X 10-° 

6.9 X 10-° 

1.5 X 10-° 

p = 10“2 

2.3 X 10“° 

2.8 X 10-'° 

2.1 X 10-° 

2.2 X 10-'° 


ber density calculation differences, but also due to differences in 
computing the stochastic GWB. 

Additionally, we computed the tension upper limits in the case 
oi p -t 1. In general, the intercommutation probability for cosmic 
superstrings can acquire values in the range p 6 [10^^, 1], depend¬ 
ing on their nature (F- or D- cosmic superstrings) ( [Jackson, Jones| 
|& Polchinski|2005| (. The effect of p 1 is the increase of the am¬ 
plitude of the stochastic GWB, without affecting the shape of the 
GW spectrum, fig* oc where k is the dependence of the scaling 
law that describes the effects of the intercommutation probability 
on the energy density of the infinite cosmic strings (poc oc 

There is no general consensus on the value of k, with differ¬ 
ent investigations suggesting a value of k = 0.6 (jAvg oustidis &| 
|Shellard|2005t or k = 1 ( |Sakellariadou|2005) . The reduced inter¬ 
commutation probability, has as a result an increased number of 
intercommutations in order to maintain the scaling evolution of the 
network, and therefore, an increased number of loops which give 
GWBs of higher amplitude. In Table we present the tension up¬ 
per limits for various cosmic superstring configurations, covering 
the whole possible range of p and k values. These upper limits, 
which are linked with small tension values, are provided by the 
networks with the smallest loop size accessible (ttcs ~ 6 x 10“"; 
see discussion in [Sanidas, Battye & Stappers|2013| (. In Fig.|I6|we 
present a set of such exclusion curves to demonstrate the change in 
the shape of the exclusion curves as we probe lower tension val¬ 
ues. The exclusion curves in the region Gp/c^ ;£ 10“'° are always 
provided by the n, = 10“* networks. The apparent discontinuities 
in some of these exclusion curves are a combinatory result of the 
abrupt local changes in the GWB upper limit curve (more evident 
in the slope region [-2,0] for the G in Fig. and our require¬ 
ment for a matching in the spectral slope of the GWB limit and the 
spectral slope of the cosmic string GW spectrum at a frequency]^ 
/= 

In Fig. [T^ we present the tension upper limits created from 
the sensitivity curve of Fig.|T^in the frequency range / e [1.7 x 
10“^, 10“^] in the case where p = 1. We did not produce results for 
the whole frequency range of Fig. since already for frequen¬ 
cies > 7 X 10“^ Hz the tension upper limits are incompatible to the 
large scale structure of the Universe as we observe it (i.e., CMB 


^ This has been verified by using smoother GWB sensitivity curves. The 
regions of the Gplc^ - Ucs parameter space where there is also a signifi¬ 
cant change in the local spectral slope of the cosmic string GW spectrum, 
might also create such discontinuities, but such were not observed in the 
various tests we have conducted to investigate this effect. If one neglects 
the requirement for a matching local spectral index, no such artifacts are 
observed. 



Figure 16. Exclusion curves for a set of network configurations with p 
With thick lines are the configurations with k = 0.6 and with thin lines the 
configurations with ^ = 1. We present exclusion curves for p = 0.1 (solid 
lines), p = 10“^ (short dashed lines) and p = 10“^ (long dashed lines). The 
dot-dashed curve is the exclusion curve for p = 1, for reference purposes. 
For all the results, we used the stochastic GWB limit of model (iii), placed 
at a frequency / = 1 yr“^. 



Figure 17. The cosmic string tension upper limits in the case of p = 1, 
obtained from the sensitivity curve presented in Fig. |13| for a range of the 
probed frequencies. For these limits, only the information about the ampli¬ 
tude of the GWB has been used. The most stringent constraint is given by 
the lowest frequency probed, whereas results for frequencies higher than 
~ 10“^, Hz have been omitted, since already for frequencies ^ 7 x 10“^ Hz 
the tension upper limits are unphysically high. 


anisotropies at large scales and galaxy distribution at small scales). 
For its computation, we used only the amplitude of the GWB per 
frequency bin, making those results less robust than those which 
also incorporate the information of the local GWB spectral index. 
The most stringent constraint is provided by the lowest frequency 

and is Gyu/c^ < l.l x 10“’. 

Our upper limit on the string tension for Nambu-Goto strings 
is slightly better than the one obtained in the extensive analysis 
performed hy the Planck Collaboration ([Planck Collaboration et al.| 


2014 1 using Planck data only {Gpjc^ < 1.5 x 10 ’), and identical 
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to that obtained when Planck data were combined with high-^ data 
from the Atacama Cosmology Telescope (ACT) and the South Pole 
Telescope (SPT). We note the significant improvement of the EPTA 
result in comparison to that presented in|Sanidas, Battye & Stappers| 
\2Q\2) , which was ~3 times less constraining than the best available 
limit at that time, set by the WMAP 7-year+ACT results. 

In terms of robustness, the CMB results are inherently more 
reliable than any GW-derived result because they depend only on 
the large scale properties of the cosmic string network (infinite 
strings), and not on the much more complicated details concern¬ 
ing the GW emission mechanism and the real cosmic string loop 
population. Our results, however, can be considered as quite reli¬ 
able upper limits, for several reasons. First, the approach we used 
to compute the GW spectra does not make any assumption about 
the fundamental cosmic string model parameters used, and is only 
subject to a few fundamental assumptions, such as the validity of 
the one-scale model and the scale-invariant evolution of the cosmic 
string network. In |Sousa & Avelino| ( |2013| l, the authors presented 
results based on a possible delay of the onset of the scaling evo¬ 
lution as the network transverses from the radiation to the matter 
era evolution. Such a scenario leads to an increase in the ampli¬ 
tude of the resulting stochastic GWB, with a consequent strength¬ 
ening of the upper limits, in comparison to our results where we 
assumed that the scaling evolution is always maintained and used 
a simple linear transition between the values of characteristic pa¬ 
rameters of the cosmic string network (i.e., the number of infinite 
strings within our horizon) from their radiation era to their matter 
era values. Therefore, our upper limits remain robust even if such a 
possibility is true. Second, we have not included the GW emission 


from kinks on the infinite strings ( Hmdmarsh|1990 [Sakellariadou 

1990 

|Allen & Shellard|1992||Kawasaki, Miyamoto & Nakayama 

2010 

1 or the emission of GWs due to the scaling evolution of the 


cosmic string network in the radiation era per se jPigueroa, F[ind-| 
|marsh & Urrestill^|2013| l. These mechanisms, even though they 
contribute to the general string GWB, can be a few orders of mag¬ 
nitude smaller than the loop emission and omitting them from our 
calculations will not alfect our upper limits. Note, however, that if 
the GWB originating from loops is not detectable by PTAs (i.e., in 
the case of very small loop creation by the network), the GWB cre¬ 
ated by the aforementioned mechanisms is the only one that can be 
detected is this frequency window. 


6.3 Relic Gravitational Waves 

Quantum fluctuations of the gravitational field in the early Uni¬ 
verse, amplified by an inflationary phase, are expected to produce a 
stochastic relic GWB (see e.g.|Grishchuk|1976||1977||Starobinsky| 
|1980[[Onde|1982) . Observations of this radiation would provide a 
unique insight into poorly understood processes in the very-early 
Universe, at energy scales ~ lO'® GeV and cosmic times ~ 10“^^ s 
( |BICEP2/Keck et al.|2015[|Ade et al.|2014^ . At long wavelengths, 
gravitational-waves generated during an inflationary epoch produce 
a characteristic signature in the polarization of the CMB radiation, 
as well as CMB temperature anisotropies ( |Grishchu^|2005| l. At 
shorter wavelengths, such as the PTA observational window, this 
radiation manifests itself as a contribution to the present day energy 
density spectrum The background spectrum is directly re¬ 

lated to the primordial tensor spectral index n, and the equation of 
state of the early-Universe (see e.g. |Zhao|201 l[[^ao et al.|2013^ . 

For standard single-field inflationary models, n, is related to 
the scalar-to-tensor ratio rhyn, = -r/&, ( [Copeland et al.|1993| l. 
If we further assume that at the end of inflation the Universe is 


not characterised by a ‘stiff’ equation of state but enters a radiation 
dominated era, w = 1/3, the spectrum of the background can be 
written as in jZhao et al.| ( |2013} : 


Ggw"(/) ~ 1-3 X 10-‘'^10‘°''' 



(58) 


where it was implicitly assumed that h = 0.6711. Current results 
from |BICEP2/Keck et al.| ( [20T5l ( set a limit of r < 0.12, and there¬ 
fore the background spectrum is almost flat over a wide frequency 
range. By comparing the frequency dependency with our model we 
can therefore see that y 5. 

We can now use the Bayesian analysis methods reported in 
Sec. |5.1.f] to perform an analysis where we vary the intrinsic noise 
parameters for the pulsars, and fix the spectral index of the corre¬ 
lated GWB term to y = 5, which implicitly assumes n, = 0, see e.g. 

. This yields a 95% upper limit of A < 1.4x10“'^ 
which translates to 


Zhao et al. 


(2013 


< 1-2 X 10-’ , 


(59) 


a factor of 9 improvement from previously reported limits in jPe 
morest et d.|p013]l, and a factor of 16 more constraining than jjenet 


et al.|(|2006^. Eq. j58|l can be inverted to yield a limit on the scalar- 


to-tensor ratio of r < 2.5 x 10® (with n, = 0). While standard in¬ 
flationary scenarios assume n, < 0, string-gas cosmology predicts 
a positive tilt in the primordial tensor spectrum (see e.g. jBranden-j 
jberger, Nayeri & Patil| ( |2014^ ; |Brandenberger et al.| ( |2007| ) which 
we cannot yet rale out. Future searches can set limits on n,, provid¬ 
ing the cosmology community independent measurements of this 
value, as well as independent constraints on r. 

The limit on < 1.2 x 10-’ is a factor ~ 10® 

from the predicted value reported in Eq. ( |58^ , which may even 
be beyond the capabilities of a future PTA using the Square 
Kilometer Array (SKA; [Janssen et al.||20f5| l, but is significantly 
more stringent than the indirect Big-Bang-Nucleosynthesis limit 
Jng„(/)d(ln/) < 1.lx 10“®(Av-3), where Av is the effective num¬ 
ber of neutrino species at the time of Big-Bang-Nucleosynthesis, 
see e.g. [Allen[ {\991\\ [Maggiore] ( [2000[ l. Current [Planck Collab-[ 


[oration et al. ( [2015) limits place a limit on A„ < 3.7. CMB 
experiments also set limits of a similar order of magnitude, 
f ngw(/)d(ln/) < 8.7x10-® 


Sendra & Smith 


(2012i; 


|Pierpaoli & KamionkowskT|([2006), whereas ground-based 


Smith, 
interfer- 

ometers, which measure the relic GWB at specific frequency inter¬ 
vals, have recently been able to do better. Indeed, LIGO and Virgo 
reported new constraints in four different frequency bands, the most 
stringent being at / = 41.5- 169.25 Hz, where - 5.6x I0“® 

at 95% confidence, with an assumed Ho - 68 km/s/Mpc, [Aasi et al.[ 
([2014). To make a direct comparison to our result, we set h = 0.68 


in Eq. 1 59 1 and find = 2.6 x 10 over three orders of mag¬ 


nitude more constraining. 

We want to stress, however, that models such as those 
described in Grishchuk j2005|; jBrandenberger, Nayeri & Patil| 


( [2014) ; [Brandenberger et al.[)2007) with values of y 6 [4.6,5] with 
n, e [0,0.9], may lead to much larger values ~ 10“'‘^, 

which may be within the reach of the SKA, see [Zhao et al.[j2013) 
for further details. 


7 CONCLUSIONS 

In this paper we have used a 6 pulsar subset of the recent EPTA data 
release presented in D15 to set a robust limit on the amplitude of a 
stochastic GWB using several models. When considering a power 
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law model for the background we obtain a limit of A = 3.0 x 10“'^ 
at a spectral index of y = 13/3, consistent with a GWB dominated 
by SMBHBs, equivalent to ng„(/)/!^ = l.lx 10“* at 2.8 nHz. When 
allowing the spectral index to vary freely over a prior range from 
0 —> 7, A = 1.3 X lO^*'^. This limit was obtained using a Bayesian 
analysis, in which we fit simultaneously for the intrinsic spin-noise 
and DM variation parameters for each pulsar, along with the GWB 
and additional common signals, including clock and Solar System 
ephemeris errors. We stress that the simultaneous fit of the GWB 
signal with the individual pulsar noise parameters, and additional 
sources of common noise is crucial to obtain a robust limit. Fix¬ 
ing the intrinsic pulsar noise to the maximum likelihood values ob¬ 
tained by the single pulsar analysis and searching for a correlated 
signal a posteriori erroneously leads to an upper limit which is a 
factor of two more stringent. A series of simulations, and a parallel 
frequentist pipeline employing the optimal statistic yields consis¬ 
tent results, corroborating the robustness of our analysis. We also 
present a more general analysis, where we do not use a power law 
model for the background, but obtain limits on the correlated power 
spectrum at a series of discrete frequencies, and show that our sen¬ 
sitivity is greatest at a frequency of ~ 5 x 10“® Hz. 

In both cases we performed model selection using the 
Bayesian evidence for models that include a common red noise 
process that is either correlated between pulsars according to the 
isotropic overlap reduction function, or that is uncorrelated be¬ 
tween the pulsars in the dataset. We obtained a difference in the 
logarithm of the Bayesian evidence of -1.0 ± 0.5 for the power law, 
and 0.2 ± 0.3 for the more general model, indicating that the dataset 
is not able to differentiate between these two cases. We confirm 
this result by obtaining confidence intervals for the correlation co¬ 
efficients between pulsars as a function of their angular separation 
on the sky and find them to be consistent both with zero correlation, 
and the Hellings-Downs curve. 

Finally, we discussed the implications of our analysis on the 
astrophysics of SMBHBs, and derived upper limits on the string 
tension of a cosmic (super)string network and for a relic GWB. Our 
upper limit of A = 3.0 x 10“'^ at a spectral index of y = 13/3 skims 
the region of the expected GWB predicted by recent astrophysi- 
cal models, but is still too high to place stringent constraints on 
the cosmological SMBHB population. An improvement of a factor 
2-3 would place our sensitivity at the heart of the expected signal 
range for the included models, placing considerable constraints on 
possible populations of merging supermassive black holes. In the 
case of a Nambu-Goto field theory cosmic string network, the upper 
limit on the string tension was evaluated to be Gpjc^ < 1.3 x 10“’, 
identical to the best so far result from CMB investigations; the re¬ 
sult presented by the Planck Collaboration combining data from 
Planck, SPT and ACT. Planck has managed to measure the temper¬ 
ature anisotropies of the CMB to an unprecedented detail, and it is 
expected that the string tension limits will not be improved signif¬ 
icantly in the future unless there is an inclusion of CMB polarisa¬ 
tion data. On the other hand, the constraints from PTAs will con¬ 
tinue to improve significantly as longer and more precise datasets 
are obtained, since p oc n'’’. Therefore, the PTA constraints on 
the string tension, as long as there is a careful treatment of all 
the involved uncertainties, have the potential to be the most strin¬ 
gent in the coming years, until full-sky CMB polarisation instru¬ 
ments become a reality (i.e., COrE-tQl. Our limit on the relic GWB, 
= 1.2 X 10“^, is a factor of 9 more constraining than 


^ https://hangai'.iasfbo.inaf.it/core/ 


the previous NANOGrav limit reported in|Demorest et al.|(|2013 


and 16 times more constraining than the last PPTA limit, see Jenet| 
|et al.| ( [2006| l. Although the expected level of the relic GW energy 
density is ~ 10“'^ in the PTA band, this number may 

increase by an order of magnitude for models with non-flat pri¬ 
mordial spectra (i.e. with nonzero tensor indices, n,), described in 
|Grishchuk| (|2005|l. Such models describe a relic GWB which may 
just be within the grasp of the SKA according to studies by |Zhao| 
[etaL] ( |2^ . 
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